arXiv:1509.01540v2 [cond-mat.mes-hall] 14 Sep 2015 


Coupling a nano-particle with isothermal fluctuating hydrodynamics: 
Coarse-graining from microscopic to mesoscopic dynamics 

Pep EspanoP, Aleksandar Donev^^’^ 

Dept. Fisica Fundamental, Universidad Nacional de Educacion a Distancia, Aptdo. 6 OI 4 I E-28080, Madrid, Spain 
Courant Institute of Mathematical Scienees, New York University 
251 Mercer Street, New York, NY 10012 
(Dated: September 15, 2015) 

We derive a coarse-grained description of the dynamics of a nanoparticle immersed in an isothermal 
simple fluid by performing a systematic coarse graining of the underlying microscopic dynamics. As 
coarse-grained or relevant variables we select the position of the nanoparticle and the total mass 
and momentum density field of the fluid, which are locally conserved slow variables because they 
are defined to include the contribution of the nanoparticle. The theory of coarse graining based on 
the Zwanzing projection operator leads us to a system of stochastic ordinary differential equations 
(SODEs) that are closed in the relevant variables. We demonstrate that our discrete coarse-grained 
equations are consistent with a Petrov-Galerkin finite-element discretization of a system of formal 
stochastic partial differential equations (SPDEs) which resemble previously-used phenomenological 
models based on fluctuating hydrodynamics. Key to this connection between our “bottom-up” and 
previous “top-down” approaches is the use of the same dual orthogonal set of linear basis functions 
familiar from finite element methods (FEM), both as a way to coarse-grain the microscopic degrees 
of freedom, and as a way to discretize the equations of fluctuating hydrodynamics. Another key 
ingredient is the use of a “linear for spiky” weak approximation which replaces microscopic “fields” 
with a linear FE interpolant inside expectation values. For the irreversible or dissipative dynamics, 
we approximate the constrained Green-Kubo expressions for the dissipation coefficients with their 
equilibrium averages. Under suitable approximations we obtain closed approximations of the coarse¬ 
grained dynamics in a manner which gives them a clear physical interpretation, and provides explicit 
microscopic expressions for all of the coefficients appearing in the closure. Our work leads to a model 
for dilute nanocolloidal suspensions that can be simulated effectively using feasibly short molecular 
dynamics simulations as input to a FEM fluctuating hydrodynamic solver. 


I. INTRODUCTION 


The study of the Brownian motion of rigid particles 
suspended in a viscous solvent is one of the oldest sub¬ 
jects in nonequilibrium statistical mechanics since the pi¬ 
oneering work of Einstein [ij. Nevertheless, it was not 
until the seventies that it was realized how subtle diffu¬ 
sion in liquids is , and to this day there remain open 
fundamental questions about the collective diffusion in 
colloidal suspensions. For example, the validity of Pick’s 
macroscopic law is questioned for suspensions confined 
to a two dimensions 0 , and it remains as a substan¬ 
tial mathematical challenge to prove that a local Fickian 
equation is the law of large numbers in three dimensions, 
even for dilute suspensions El- These questions are not 
of purely academic interest since diffusion is of crucial 
importance in a number of applications in chemical en¬ 
gineering and materials science, such as the study of the 
dynamics of passive or active |l2l . Il3| particles in suspen¬ 
sion, the dynamics of biomolecules in solution |l4, Il5l | , 
the design of novel nanocolloidal suspensions |l6l4l8j , and 
others. The importance of coarse-graining to the study 
of diffusion in nanocolloidal suspensions is easy to ap¬ 
preciate; the number of degrees of freedom necessary to 
simulate Brownian motion directly using Molecular Dy¬ 
namics (MD) is large enough to make this approach pro¬ 
hibitively expensive. In this paper, we derive from “first 
principles” a coarse-grained dynamic equation for the po¬ 


sition of a nanoparticle immersed in a simple fluid, fully 
taking into account hydrodynamic effects. 

The key source of difficulty in the theoretical and com¬ 
putational modeling of colloidal diffusion is the presence 
of viscous dissipation in the surrounding fluid. This 
hydrodynamic dissipation in the solvent induces long- 
ranged hydrodynamic fields that couple the motion of 
the solute particles to boundaries and to other parti¬ 
cles. These effects are termed hydrodynamic interac¬ 
tions in the literature, but it should be kept in mind 
that these “interactions” are different in nature from di¬ 
rect interactions such as steric repulsion or long-ranged 
attractions among the colloids. The well-known Smolu- 
chowski or Brownian Dynamics (BD) [l^ 10 approach 
captures the effect of the solvent through a mobility ma¬ 
trix that is approximated using hydrodynamic models 
based on assumptions that are of questionable validity 
for nanoscopic particles. In particular, a gold nanocol¬ 
loid and a biomolecule such as a protein can only be 
distinguished in BD based on an effective hydrodynamic 
no-slip surface but not based on the nature of their inter¬ 
action with the solvent. This makes BD unsuitable for 
capturing multiscale effects such as slip on the surface 
of the particle, layering of the solvent molecules around 
the colloid, transient hydrogen bond networks around the 
protein, etc. 

The fluctuation-dissipation balance principle informs 
us that viscous dissipation is intimately related to fluctu- 
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ations of the fluid velocity. It is well-known that diffusion 
in liquids is strongly affected by advection by thermal ve¬ 
locity fluctuations i, and that nonequilibrium 

diffusive mixing is accompanied by “giant” long-range 
correlated thermal fluctuations |23-l^. As explained in 
detail in Refs. (ill. [28l - l^ . there is a direct relation be¬ 
tween these unusual properties of thermal fluctuations in 
liquid solutions and Brownian Dynamics. Specifically, a 
simplified model of colloidal diffusion based on incom¬ 
pressible fluctuating hydrodynamics can be mapped one- 
to-one to the equations of BD and related Dynamic Den- 
sity Functional Theories (DDFT) with hydrodynamics 
|31L 1321 ; this deriyation shows that hydrodynamic inter¬ 
actions are nothing more nor less than hydrodynamic cor¬ 
relations induced by the thermal fluctuations in the sol- 
yent. Such a fluctuating hydrodynamic model (Til. [28l - [30j 
explains the appearance of giant nonequilibrium fluctua¬ 
tions in the concentration of colloidal particles, justifies 
the Stokes-Einstein relation in the limit of large Schmidt 
numbers [sj, and describes the important influence of 
boundaries in confined suspensions |23l l29j . If one wants 
to further account for inertial effects and compressibility 
of the fluid, as crucial for modeling the effect of ultra¬ 
sound on colloidal particles [sJl or the acoustic vibrations 

S duced by suspended particles (^ or micro-organisms 
, one can use a similar model but describe the fluid us¬ 
ing compressible fluctuating hydrodynamics (^. l3^ . 

In this work we consider coupling compressible isother¬ 
mal fluctuating hydrodynamics to a suspended nanocol- 
loid al p article. Unlike previous phenomenological mod¬ 
els |llL [^. we obtain our equations from 

the underlying microscopic dynamics by using the The¬ 
ory of Coarse-Graining (TCG) as developed by Green 
[43 and Zwanzig (also see the textbook (i^), to¬ 
gether with a sequence of careful approximations that 
preserye the correct structure of the exact (but formal) 
coarse-grained equations. Our deriyation is important 
for seyeral reasons. Firstly, our work proyides a micro¬ 
scopic foundation for the types of models used in exis ting 
theoretical and computational work [llL 1^ . [ 3 ^ 

m. Secondly, and more importantly, our deriyation 
leads to microscopic Green-Kubo type formulas for the 
transport coefficients that appear in the coarse-grained 
equations. This allows for these coefficients to be esti¬ 
mated from molecular dynamics computations, thus fully 
taking into account microscopic effects that are difficult 
if not impossible to include in purely continuum mod¬ 
els. Thirdly, our deriyation will lead us to first con¬ 
struct a microscopically-justified fully discrete form of 
compressible isothermal fluctuating hydrodynamics that 
is second-order accurate while also maintaining discrete 
fluctuation-dissipation balance to second order. 

This last contribution is in itself a significant exten¬ 
sion of prior work [i^, fully consistent with the ap¬ 
proach to nonlinear fluctuating hydrodynamics proposed 
in our recent work (^ . Specifically, the coarse-grained 
equations we deriye here by following a “bottom-up” ap¬ 
proach can also be deriyed by a “top-down” approach in 


which one starts from a (phenomenological) system of 
formal stochastic partial differential equations and ap¬ 
plies a Petrov-Galerkin hnite-element discretization . 
Our work therefore proyides a direct and explicit link 
between the microscopic discrete dynamics and meso¬ 
scopic continuum fluctuating hydrodynamics. The phys¬ 
ical insight that is necessary to construct phenomeno¬ 
logical fluctuating hydrodynamics equations translates 
in this paper into physical insight required when con¬ 
structing suitable approximations or closures of a num¬ 
ber of intractable microscopic expressions. The “bottom- 
up” procedure clearly reveals all of the required terms 
in the coarse-grained equations and provides microscopic 
expressions for the required coefficients. 

At first sight, it may seem like the equations of Smolu- 
chowski that underlie Brownian dynamics have a well- 
known microscopic derivation. Indeed, it is not difficult 
to construct a text-book TGG for the dynamic equa¬ 
tion describing the positions of the colloidal particles 
( 4 ^ . This leads to the well-known expression for the 
hydrodynamic mobility (diffusion tensor) as the time in¬ 
tegral of the correlation function of the velocities of the 
solute particles, conditional on the particle’s positions. 
It should, however, quickly be recognized that this well- 
known expression, while correct, is not useful in prac¬ 
tice, for several reasons. Firstly, this integral must be 
computed anew for every configuration of the suspended 
particles. Secondly, even if one could run a new MD 
calculation at every step in a BD simulation, it is im¬ 
portant to realize that these MD computations are un¬ 
feasible in practice because they must be very long on 
microscopic scales. Namely, it is well-known that the 
slow viscous (diffusive) dissipation of momentum in the 
fluid makes the velocity correlation functions have long 
(power-law) hydrodynamic tails; it is the integral of these 
tails that gives the hydrodynamic correlations (interac¬ 
tions) among the particles, as well as finite-size effects on 
the diffusion coefficient for confined particles [IJ. There¬ 
fore, to correctly capture hydrodynamic effects the time 
integral in the Green-Kubo expression for the diffusion 
tensor must extend to at least the time it takes for mo¬ 
mentum to diffuse throughout the whole system; while 
this time is typically short compared to the time scale at 
which the solute particles move, it is very long based on 
MD standards. 

By contrast, in the equations derived here the Green- 
Kubo integrals can be computed via feasible (short) MD 
simulations. This is because all of the hydrodynamics, 
such as the effects of sound [s^ or viscous dissipation 
( 3 ^ are captured by explicitly resolving the (fluctuating) 
hydrodynamics of the solvent using a grid of hydrody¬ 
namic cells, and only the remaining local and short-time 
effects need to be captured by the microscopic simula¬ 
tions. In the present work, we consider suspensions that 
are sufficiently dilute to allow us to neglect the direct 
(as opposed to hydrodynamic) interactions among the 
colloids and focus our derivation on a single particle im¬ 
mersed in a viscous liquid; hydrodynamic interactions 
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among the particles are still captured because they are 
mediated by the explicitly resolved surrounding fluid dy¬ 
namics. In fact, we believe that in many cases of interest 
the coarse-grained diffusive dynamics can effectively be 
simulated by a priori performing a small number of short 
MD simulations of a single particle in a small (say peri¬ 
odic) domain. Crucial to the above is the fact that in the 
present work the hydrodynamic cells are assumed to be 
significantly larger than the nanoparticle itself. 

In the next section we explain in more detail the basic 
assumptions and thus limitations of our model. Briefly, 
our model assumes that the solvent is a simple isotropic 
single-component fluid. We do not explicitly consider 
energy transport and thus limit our work to isothermal 
suspensions. We only consider dilute suspensions of nano 
particles. The extension to denser suspension leads to a 
significantly more complicated theory of liquid mixtures 
that is well beyond the scope of this work. The lim¬ 
itation to nanoscopic particles is not essential and the 
equations developed here can be used also for larger par¬ 
ticles such as micron-sized colloids; however, in this case 
the MD simulations required to obtain the values of the 
Green-Kubo integrals that appear in the coarse-grained 
equations would again become unfeasible and a different 
approach is advised. We will also assume that the parti¬ 
cle is effectively spherical so that describing the position 
of its center of mass is sufficient without requiring us to 
also resolve its orientation. Our theory assumes a separa¬ 
tion of time scales between the positions of the particles 
and their velocities, and we do not include the velocities 
of the colloidal particles in the description. More pre¬ 
cisely, it requires that the Schmidt number of the solute 
particles be very large. This is not a significant limita¬ 
tion in practice since the Schmidt number of even a sin¬ 
gle solvent molecule is typically very large in liquids. In 
particular, our theory can be used to describe collective 
diffusion of tagged solvent particles (i.e., self-diffusion). 

In Section m we explain the basic notation and con¬ 
cepts, and carefully select and define the coarse-grained 
(slow) variables in terms of the microscopic degrees of 
freedom. We then proceed to carefully examine the re¬ 
versible (non-dissipative) part of the dynamics. In par¬ 
ticular, in Section BlIAl we give exact results that are not 
useful on their own right since they lead to equations that 
are not closed explicitly . However, by making a series of 
approximations based on a key “linear for spiky” approx¬ 
imation we are able to derive an approximate closure for 
the reversible dynamics in Section riIIBl In Section HVl we 
apply the same approximation to the irreversible (dissi¬ 
pative) part of the dynamics, together with another im¬ 
portant approximation in which we replace constrained 
Green-Kubo expressions with unconstrained equilibrium 
Green-Kubo averages. The key results of our calcula¬ 
tions are then collected and discussed in Section El We 
first give an approximate but closed form for the coarse¬ 
grained discrete dynamics, and then discuss the relation 
of these discrete equations to continuum models in Sec¬ 
tion IV Gl A comparison of our results to phenomenologi¬ 


cal models and a discussion of their significance and range 
of validity is given in Section IVII A number of technical 
calculations are detailed in an extensive Appendix. 


II. COARSE-GRAINING 


In this section, we give the basic ingredients required to 
perform the coarse-graining of the microscopic dynamics 
for our specific system. We begin with a general overview 
of the theory and then specialize to the case of a nanopar¬ 
ticle suspended in a simple liquid by explaining the de¬ 
tails of the microscopic dynamics and the definition of 
the coarse-grained variables. 


A. The Theory of Coarse-Graining 


In this section, we review the theory of Goarse- 
Graining or Non-Equilibrium Statistical Mechanics as es¬ 
tablished by Green and Zwanzig [ij . The theory al¬ 
lows to construct the dynamic equations for the probabil¬ 
ity distribution of a set of coarse-grained (GG) variables 
that describe the state of a system at a coarse level of de¬ 
scription. The theory states that, under the assumption 
that the CG variables are sufficiently slow as compared 
with the eliminated degrees of freedom, the system fol¬ 
lows a diffusion process in the space of GG variables. 
The resulting dynamic equation for the probability dis¬ 
tribution of the GG variables is given by a Fokker-Planck 
equation (FPE), where both the drift and diffusion terms 
are given in microscopic terms. 

The coarse-grained variables are selected functions 
x{z) in phase space, i.e. they depend on the set of posi¬ 
tion and momenta z of the molecules of the system. We 
follow the convention that a hatted symbol like x{z) de¬ 
notes a function in phase space that may take numerical 
values X. The selection of the relevant variables x{z) is a 
crucial step in the description of a non-equilibrium sys¬ 
tem. A crucial requirement is that they are slow variables 
( 4 ^ . When this is the case, the probability distribution 
of a set of relevant variables x obeys the FPE 


dtP{x, t) 


A 

dx 


f)nj 

A{x)-V{x)~{x) 


P{xA) 


+ ksT 


dx 


'P{x)~P{x,f) 


( 1 ) 


The different objects in this equation have a well-defined 
microscopic definition. For example, the reversible drift 

is 


^(x) = {LxY (2) 


where L is the Liouville operator and the conditional ex- 
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pectation is defined by 

^ • ■)" = / dzp<^^{z)5{x{z) -x)--- (3) 

where p’^'^{z) stands for the microscopic equilibrium dis¬ 
tribution and 6(x(z) — x) is actually a product of Dirac 
delta functions, one for every function x{z). The equilib¬ 
rium distribution of the relevant variables is 

P^^{x) = f dzp^^[z)6{x{z) — x) (4) 


and is closely related to the bare free energy of the level 
of description x which is defined through 

n{x) = -kBT\nP^^{x) (5) 


Here kg is Boltzmann’s constant and T the temperature 
of the equilibrium state. We will refer in this work to the 
bare free energy also as the coarse-grained Hamiltonian 
because of the particular form that T-L^x) acquires at the 
hydrodynamic level of description. When non-isothermal 
situations are considered one rather introduces the en¬ 
tropy of the level of description as S{x) = fcs lnP®'^(a;), 
according to Einstein formula for fluctuations. 

Finally, the symmetric and positive semidefinite (45| 
dissipative matrix ^{x) is the matrix of transport coeffi¬ 
cients expressed in the form of Green-Kubo formulas, 


V{x) = 


keT . 


{QLX exp{iQLt'}QLX)^ dt' (6) 


The term QLx is the so called projected current. The 
projection operator Q is defined from its action on any 
phase function B{z) 


QB{z) = B{z) - (7) 


The dynamic operator exp{iQLt'} is usually named the 
projected dynamics, which is, strictly speaking different 
from the real Hamiltonian dynamics exp{Lt'}. The pro¬ 
jected dynamics can be usually approximated by the real 
dynamics but, in order to avoid the so called plateau prob¬ 
lem [ 4 ^, then the upper infinite limit of integration in 
Eq. (ini) has to be replaced by r, a time which is long in 
front of the correlation time of the integrand, but short 
in front of the time scale of evolution of the macroscopic 
variables (dR IdOl - f^ , this is 

P(x) = -^/ {QLXexp{iLt'}QLX)^dt' (8) 
ksT Jq 

In general, it is expected that different elements of the 
matrix may require different values of r. 

The Ito stochastic differential equation (SDE) that is 
mathematically equivalent to the FPE ([1]) is given by 

§ = Mx) - + kBT^ Wi] + §(i) (9) 


where ^(x) = B{x) is a linear combination of white 
noises, formally time derivatives of a collection of in¬ 
dependent Wiener processes (Brownian motions) B{t), 
where the amplitudes satisfy the Fluctuation-Dissipation 
Balance (FDB) condition 

B{xfB{x) = 2kBTV{x) (10) 

In summary, the three basic objects that determine the 
dynamics (either in the FPE ([T|) or the SDE ([H]) forms) 
and that need to be computed in the theory are the bare 
free energy %{x), the reversible drift .4(x), and the dis¬ 
sipative matrix P(x). 

The reversible drift can also be written in the form (45| 

A^{x) = Lf^„{x)-^{x) - kBT (x) (11) 

where the skew-symmetric reversible matrix is defined as 

L^4x) = ({X^,X,}f (12) 

where {•,•} is the Poisson bracket. Here and in what 
follows, Einstein convention that sums over repeated in¬ 
dices is assumed. Note that the form of the drift m 
ensures automatically the Gibbs-Boltzmann distribution 
P®'^(a;) cx is the equilibrium solution of (H)), even 

for approximate forms of the reversible matrix L(x) and 
the CG Hamiltonian TL{x), and, thus, is the preferred 
form for the reversible drift in the present work. 


B. Selection of Coarse-Grained Variables 

The most important step in the TGG is the selection 
of the relevant (coarse-grained) variables. This selection 
must be guided by physical intuition and the presence 
or absence of separation of time scales. The key guid¬ 
ing principle is that the relevant variables must evolve 
much more slowly than all other variables that cannot 
be expressed entirely in terms of the relevant variables. 
This allows us to make a Markovian approximation of 
the coarse-grained dynamics, which takes the form of a 
Fokker-Planck equation for the probability distribution 
of relevant variables, or equivalently, of a stochastic dif¬ 
ferential equation for the instantaneous (fluctuating) rel¬ 
evant variables. 

Ultimately, one is often only interested in the positions 
(and possibly orientations) of the colloidal particles, elim¬ 
inating the solvent from consideration entirely. This is 
possible to do via TGG because indeed in liquids mass 
diffusion is very slow compared to momentum and heat 
diffusion, and thus the positions of the particles are much 
slower than the hydrodynamic fields. Indeed, following 
the TGG using only the positions of the particles leads 
to the well known equations of Smoluchowski or Brown¬ 
ian dynamics, with well-known Green-Kubo expressions 
for the hydrodynamic mobility (equivalently, diffusion) 
matrix (see, for example. Section V in [i^). As we ex- 







5 


plained above, this level of description is not sufficiently 
detailed to allow us to describe a number of important 
microscopic effects that occur in the vicinity of the parti¬ 
cle surface. While in the present work we do not capture 
explicitly the slip at the surface and the layering effects 
around a nanoparticle, we do take into account such ef¬ 
fects implicitly through the microscopic expressions that 
enter in the theory. Furthermore, the Green-Kubo for¬ 
mulas for the mobility are not useful in practice and one 
must close the equations by using a pairwise approxima¬ 
tion to the mobility matrix based on far-held expansions 
for Stokes flow. 

To go to a more fundamental (microscopically more 
informed) level of description we must include solvent 
degrees of freedom as well. We want to describe the sol¬ 
vent molecules at the hydrodynamic rather than the mi¬ 
croscopic level since it is not reasonable to keep track 
of the positions and momenta of every molecule in the 
system. At macroscopic scales, a fluid appears as a con¬ 
tinuum that is described with smooth helds obeying the 
well-known Navier-Stokes equations. The “held” concept 
is tricky, though, because a held is a mathematical ob¬ 
ject that has inhnitely many degrees of freedom, while 
the actual huid system has a hnite number of degrees of 
freedom. Of course, the helds are dehned above a cer¬ 
tain spatial resolution much larger than the typical size 
and distances between molecules of the huid. At these 
macroscopic scales the held at one point of space effec¬ 
tively represents a very large number of molecules that 
move in a coherent manner. When one descends down 
to mesoscopic scales, molecules do not move that coher¬ 
ently, and one starts appreciating the discrete nature of 
the huid. In other words, the average behavior and the 
actual behavior of the huid molecules start to differ, and 
it is necessary to describe a huid system with hydrody¬ 
namic equations that are intrinsically stochastic. The 
hrst phenomenological theory for such huctuating hydro¬ 
dynamics was proposed by Landau and Lifshitz, who in¬ 
troduced the concepts of random stress and heat huxes, 
to be added to the usual Newtonian stress and Fourier 
heat hux [s^. 

From a mathematical point of view, the nonlinear 
stochastic partial diherential equations (SPDEs) of huc¬ 
tuating hydrodynamics are ill-dehned. In other words, 
a continuum limit of sequences of more rehned other¬ 
wise reasonable discrete versions of the partial differen¬ 
tial equation does not exist. From a physical point of 
view, though, this is not much of a problem because 
we know that the continuum limit cannot be realized 
without hrst encountering the atomistic nature of mat¬ 
ter. For these reasons, it is necessary to dehne discrete 
hydrodynamic variables by averaging over a number of 
nearby molecules, and use these discrete variables in the 
TCG. In this work, following the approach developed in 
a sequence of prior works (^ . 1 ^ 1^ , we dehne discrete 
hydrodynamic helds by placing a hxed (Eulerian) grid 
of hydrodynamic nodes and associating to each node a 
huid density and momentum averaged over a hydrody¬ 


namic cell associated to that node. In the present work 
we compute with more rigor some of the conditional ex¬ 
pectations that were plausibly approximated in (d^ . In 
order to have a reasonable hydrodynamics description we 
need to have hydrodynamic cells that contain many sol¬ 
vent molecules; here we consider simple liquids for which 
hydrodynamic cells containing many molecules will also 
be much larger than the mean free path. 

For a colloidal particle that is much larger than the 
solvent molecules, the hydrodynamic how around the 
nanoparticle can be resolved with small (compared to the 
size of the nanoparticle) hydrodynamic cells that, nev¬ 
ertheless, still contain many solvent molecules. In this 
situation, the discrete fluid mass density and the dis¬ 
crete fluid momentum densities g^, where p indexes the 
hydrodynamic nodes, would only include contributions 
from the solvent particles. At such a level of description 
it is necessary to include both the position R and the 
momentum P of the nanoparticle in the list of relevant 
variables because even though P is much faster than the 
position, it evolves on the same time scale as the hydro- 
dynamic momentum around the particle. This level of 
description has been traditionally used for the descrip¬ 
tion of Brownian motion of colloidal particles coupled 
with fluctuating hydrodynamics ii. We do not con¬ 
sider this case here; for a phenomeno log ical model of this 
type we refer the reader to Refs. [^[^1^. It is impor¬ 
tant to note that it is inconsistent to keep the velocities 
and thus inertial dynamics of the particles without also 
accounting for the viscosity and inertia of the surround¬ 
ing fluid. This is because there is not a separation of 
time scales between the velocities of the particles and 
the velocity of the surrounding fluid; the only consistent 
coarse-grained implicit-fluid level of description is that of 
Brownian dynamics, as explained in detail by Roux j^. 

Here we consider a nanoparticle that is not much 
larger than the fluid molecules, so that the hydrody¬ 
namic cells are much larger than the nanocolloidal par¬ 
ticle, i.e., we have a “subgrid” colloidal particle. In par¬ 
ticular, the “nanoparticle” particle could be just a tagged 
fluid molecule when modeling self-diffusion in a liquid. 
Since the momentum of the particle evolves on the same 
time scale as the solvent molecules with which it collides, 
more precisely, since the fluctuations of the relative ve¬ 
locity of the colloid are fast compared to hydrodynamic 
time scales, we define the hydrodynamic mass and mo¬ 
mentum density fields to include the nanoparticle con¬ 
tribution. In summary, the level of description that we 
consider in this work is characterized by the position of 
the colloid R, the (total, i.e., including the contribution 
from the nanoparticle) discrete mass density and the 
(total) discrete momentum density g^, where p indexes 
the hydrodynamic nodes. 

We make use of the standard TGG of Zwanzig where 
all the terms (CG free energy, drift, and diffusion matrix) 
are given in microscopic terms (idl . . This allows one 
to obtain the general structure of the dynamics. However 
in order to find tractable results it is crucial to make a 
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number of assumptions. All the approximations that we 
consider rely on the fact that the cells used to define the 
hydrodynamic variables are much larger than the typi¬ 
cal intermolecular distances in such a way that every cell 
contains many molecules of the fluid. In particular, we 
assume that the microscopic local density field which is 
of the form miS{r — c[i) gives, once inside conditional 

expectations, the same result as the interpolated discrete 
density variables (see Eq. (H51) below and Fig. [3]). This 
is only plausible if, again, there are many molecules per 
cell and the values of the discrete variables in neighbor¬ 
ing cells are very similar. While this is statement about 
the flow regimes for which the resulting equations apply, 
it is also an statement about the size of the fluctuation of 
the hydrodynamic variables. They need to be small, oth¬ 
erwise, the value in neighbor cells could be very different 
just by chance. In other words, the number of molecules 
per cell must be sufficiently large in order for the relative 
fluctuations to be sufficiently small. In the end, the va¬ 
lidity of the approximations made and the utility of the 
final equations we obtain can only be judged by a com¬ 
putational comparison to the true microscopic dynamics 
(molecular dynamics). 


C. Microscopic Dynamics 

In the present work we consider a simple liquid system 
of -I- I particles described with the position and mo¬ 
menta of their center of mass (see Fig. [T]for a schematic 
representation), in a periodic box. We distinguish parti¬ 
cle i = 0 as the nanoparticle which has a mass mg, typi¬ 
cally larger than the mass to of a solvent particle. At the 
microscopic level the system is described by the set z of 
all positions and momenta pi = TO^Vi (i = 0,1, • • • , N) 
of the particles. The microstate of the system evolves ac¬ 
cording to Hamilton’s equations with Hamiltonian given 
by 


H{z) 

U{q) 


W°\q) 


Po 

2 too 2m 

i=l 

N 


N 2 

i: It+'>(■<) 


+ *“(».)+ 1>"(qo) 


2=1 


N 




(13) 


i,j=l 


We have assumed a pairwise potential energy (f'iQij) 
between liquid molecules i,j separated a distance qij. 
i/®°*(g) is the potential energy of the solvent in the ab¬ 
sence of the nanoparticle, <I>™*(g) is the potential of in¬ 
teraction of the i-th solvent particle with a nanoparti¬ 
cle a distance q away, and <h®’’'*(qo) is an external time- 
independent potential acting on the nanoparticle. The 
system is assumed to have periodic boundary conditions. 

Under the assumption that the Hamiltonian is mixing. 



FIG. 1: Schematic representation of a nanoparticle (in brown) 
surrounded by molecules of a simple liquid solvent (in blue). 
Also shown is the triangulation that allows to define the dis¬ 
crete hydrodynamic variables at the nodes (in red). The 
shaded area around node jj, located at is the support of the 
finite element function '(/’^(r) and defines the hydrodynamic 
cell. 


the dynamics will sample at long times the molecular 
ensemble (S^ l given by 

(14) 

where Pq and Eq are the initial total momentum and en¬ 
ergy of the system. We will assume that in the thermody¬ 
namic limit the molecular ensemble can be approximated 
by the canonical ensemble 

peq(^) = ^exp{-/3IT(z)}, (15) 

where j3 = l/^ksT), and we use the canonical ensemble 
in the theory for simplicity. 


D. Definition of Coarse-Grained Variables 


The first step in the Theory of Coarse-Graining is to 
specify the relevant variables in terms of the microscopic 
state z of the system. In the present case, we choose as 
relevant variables the position of the nanoparticle 


R(z) = qo, (16) 

and the mass and momentum hydrodynamic “fields”. As 
we will consider fluctuations in the hydrodynamic vari¬ 
ables, the latter need to be defined in discrete terms (d^ . 
This is, we want to look at the mass and momentum 
of collections of molecules that are in a given region of 
space. To this end, we seed physical space with a set of 
M nodes, located at the points r^. Usually, the nodes are 
arranged in a regular lattice, but this is not necessary in 
what follows and arbitrary simplicial grids can be used 
(see Fig. [T]for a schematic representation). 

We define the mass and momentum densities of the 
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node /i according to 


N 

i^O 

N 

kt^{z) (17) 

i=0 

where the index i = 0 labels the nanoparticle. The basis 
function Sfj,(r) is a function (with dimensions of inverse 
of a volume) that is appreciably different from zero only 
in the vicinity of r^. This region is referred to as the 
hydrodynamic cell of node fi. We may regard the ba¬ 
sis function 6^{r) as a “discrete Dirac delta function”. 
Its specific form is discussed below. Note that both the 
mass and momentum densities contain the nanoparticle 
in their definition. It is convenient to introduce also the 
hydrodynamic fields of the solvent 

N 

i=l 

N 

= (18) 

i=l 

that do not contain in its definition the contribution of 
the nanoparticle (i.e. the particle i = 0 is excluded in the 
sum). 

We may express the discrete hydrodynamic variables 
(inj and (fT5|) in terms of the usual microscopic densities 


N 

Pv{z) = ^TOi(5(r - qi), 
N 

Sr(z) = ^Pi(5(r-qi), 
2=0 


as simple space integrals, 

Pm(z) = J drSf,{r)priz), 
Stj.{z) = J drSi,{r)gr{z), 


N 

pT\z) = '^m,S{r-q_i) 
2=1 
N 

gr'(^) = Ep.5(r-q,) 

2=1 

(19) 

P’i^\z) = J drSf,{r)pl°\z) 

if{z) = j dv5,{v)^r\z) 

( 20 ) 


Note that the two sets of variables {R, p, g} and 
{R, ^soi gsoij expressible in terms of each other. 

While we have that the densities are related as 

pf[z) = Mz) - mo^^R) (21) 

there is no way to express the momentum g as a function 
of R,g®°*. Therefore, the dynamic equations to be 
obtained for each set of variables are essentially different 
and cannot be obtained from each other through a simple 
change of variables. In other words, the two sets of rel¬ 



FIG. 2: The finite element basis function '(pij.ijc) in two dimen¬ 
sions. 


evant variables lead to physically different descriptions. 
Since the slowness of the hydrodynamic variables arises 
from the underlying conservation laws, and only the total 
mass and momentum fields are conserved quantities, the 
appropriate variables for the TCG are our chosen vari¬ 
ables {R, p, g}. 


E. The basis functions 

The actual form of the discrete Dirac delta function 
i5^(r) needs to be specified. One possibility is to use the 
characteristic function (divided by the volume of the cell) 
of the Voronoi cell of node p. For pffz) this will give the 
total mass (per unit volume) of the particles that happen 
to be within the Voronoi cell p. As we discussed in Ref. 
[^ . though, this selection is unsuited for the derivation 
of the equations governing discrete hydrodynamics from 
the Theory of Coarse-Graining. This is because the gra¬ 
dient of the characteristic function of the Voronoi cell is 
singular and leads to ill-defined Green-Kubo expressions. 
It was suggested to instead use the Delaunay triangula¬ 
tion associated with the set of nodes as a grid of finite 
elements (FE), and take the discrete delta function to be 
the linear FE basis function associated with node 

p, which has the characteristic shape of a tent in one di¬ 
mension, a pyramid in two dimensions (as shown in Fig. 
ED, and more generally a (d -|- I)-dimensional simplex in 
d dimensions. Note that the use of a Voronoi/Delaunay 
tessellation is not required, and any simplicial grid (i.e., a 
triangular grid in two dimensions or a tetrahedral grid in 
three dimensions) whose vertices are the set of hydrody¬ 
namic nodes can be used equally well (but for numerical 
purposes the grid should be kept as close to uniform as 
possible).! 

In recent work (^ . , we have argued that an even 

better selection (in terms of numerical accuracy) is given 
by a basis function S^(r) that is a linear combination of 
the (dimensionless) finite element linear basis functions 
functions '!/i’^(r) 


S^{r) = (22) 

The crucial requirement is that these basis functions are 



mutually orthogonal 


(23) 

where we have introduced double bars to denote integra¬ 
tion over space, this is 

1/1 = j d^fir) (24) 

for an arbitrary function /(r). Note that from (I22II and 
(l23)) it follows the explicit matrix form 

( 25 ) 

If we introduce the usual “mass matrix” of the finite ele¬ 
ment method 

Mtu = UM ( 26 ) 

the orthogonality condition implies that in (1^^ is 
given by the inverse of this is 

MtMa = (27) 

The basis function (r) may be regarded as a way of 

discretizing a field a(r) according to = li5^a||. The 
basis function 'i/'m(’') permits to construct interpolated 
fields out of the discrete fields a(r) = The 

orthogonality condition (1^51) ensures that if we discretize 
an interpolated field, we recover the original discrete val¬ 
ues, i.e. Ii5;^a|| = a^. This is the main motivation to use 
the slightly more involved basis function i5;_j(r) instead 
of the finite element for the definition of the CG 

variables. It turns out that this complication pays off, as 
the resulting hnite difference operators are second order 
accurate approximations of the correspondin g co ntinuum 
differential operator, even in irregular grids |56l |. 

The finite element linear basis functions satisfy a par¬ 
tition of unity and give linear consistency, 

^V'A*(r) = l, ^r^^/>^(r)=r (28) 

As a consequence of these properties, the conjugate basis 
functions i5^(r) satisfy 

XI = 1, X V,xr;.(5^(r) = r (29) 

fl fl 

where is the volume of the hydrodynamic cell 

J dr'ipf,{r) (30) 

Note that we have 

J drSf,{r) = 1, J dr r(5^(r) = (31) 


as can be proved by using (1^51) and the orthogonality 
([23]). These properties justify to call <i;j(r) a discrete 
Dirac delta function. 

The partition of unity reflected in (1^^ implies 

Xv^V(5^(r)=0 (32) 

which we will use often in proving that the resulting dy¬ 
namic equations are conservative. In fact, we define the 
total mass and total momentum of the system at the CG 
level through, 

Mt = X = X 

H i 

Pt = '^Vf,gf,{z) = Y^p, (33) 

fl i 

which are, indeed, the total mass and momentum. These 
quantities are conserved by the microscopic dynamics 
and need to be conserved by the coarse-grained dynam¬ 
ics. 

It is convenient to introduce also the following regular¬ 
ized Dirac delta function 

A(r, r') = ,5^(r)V'^(r') = A(r', r), (34) 

which is closely related to what is called the discrete 
Delta function or interpolation kernel in [2^, I37l - l4ll | . 
This function is different from zero only for distances of 
the order of the size of the hydrodynamic cells. In the 
limit of zero lattice spacing A(r,r') converges in weak 
sense to (5(r — r'). Therefore, A(r, r') can be understood 
as a Dirac delta function regularized on the scale of the 
grid. 

The regularized Dirac delta satisfies the exact identi¬ 
ties 

y dr'A(r,r')(5^(r') = (5;i(r) 

J dr'A{r,r')tjj^{r') = tp^{r) (35) 

One of the basic approximations that we will make in the 
present work is the smoothness approximation 

J dr'A(r')A(r', r) = ||A(5^||^/>^(r) ~ A(r) (36) 

for a smooth function A(r). For smooth functions the 
regularized Dirac delta acts like a Dirac delta. The ap¬ 
proximation (I36p is an exact identity for linear functions 
A(r) = a -I- r-b. Therefore, the errors committed when 
using the approximation (1361) for smooth functions are of 
second order in the lattice spacing. Sometimes, we will 
use the above identity in the form 

lA^^II U,B\\ ez \\AB\\ 


(37) 
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for any two smooth functions A(r),_B(r). 

Finally, note that one property that is not satisfied by 
the regularized Dirac delta function, as opposed to the 
Dirac delta is the following symmetry 

= ( 38 ) 

If the regularized Dirac delta function was translation- 

ally invariant, i.e. A(r,r') = A(r — r'), this would be 
obviously true. In this case, we would have in addition 
to (|35l) also the following relations, 

j dr'A{r,r')VSf,{r') = VS^{r) 

J dr'A{r, r') VV/i(r') = (39) 

Even though these identities are not fulfilled, we will as¬ 
sume that they are reasonable approximations, particu¬ 
larly if both sides are multiplied with “smooth discrete 
fields”, i.e. 

dr'A(r, r')V'a(r') ~ Va(r) (40) 

For a sufficiently smooth field a(r), the length scale of 

variation of Va(r) is much larger than the length scale 
of variation of A(r,r') and, therefore, A(r, r') acts as an 
ordinary Dirac delta. 

F. Notation 

The notation in the present work is unavoidably dense 
because many different mathematical objects need to be 
carefully distinguished. Below we present a summary of 
the notation for the case of the mass density variable 
alone. Similar symbols are used for the velocity and mo¬ 
mentum density variables. In general, hatted symbol like 
in 

N N 

PtJ.(.z) = Pr{z) = ^miS{qi-r), (41) 

2=0 2=0 


denote phase functions. The numerical values taken by 
a phase function are denoted without hat as in, for ex¬ 
ample, p^. The subscript is used here to distinguish the 
specific node p for discrete variables such as or the 
specific point in space for continuum fields such as pr. 
Overlined symbols like 

p(r) = ^/>^(r)p^ (42) 

denote continuum fields which are interpolated from dis¬ 
crete “fields”. Differential operators act only on the sym¬ 
bol immediately to their left unless otherwise indicated 
by parenthesis, dot denotes contraction, and colon a dou¬ 
ble contraction. 

III. THE REVERSIBLE DRIFT 

In this section, we present a number of exact and then 
approximate results for the reversible part A{x) of the 
dynamics and the bare free energy 'H{x) for the present 
level of description. 

The exact results presented in section IIII Al are ob¬ 
tained by integrating the microscopic momenta in the 
microscopic definitions ((2) and (|4]) for these quantities. 
This integration is possible because we assume that the 
equilibrium ensemble is given by the canonical ensemble 
(fT51l and the resulting space integrals involve relatively 
simple Gaussian integrals of the kind discussed in Ap¬ 
pendix [E] The molecular ensemble (03 can also be used 
at the expense of much cumbersome expressions. We as¬ 
sume that in the thermodynamic limit both ensembles 
are equivalent and we opt for the simpler case. In Sec¬ 
tion IIII Bl we approximate the exact results in order to 
obtain a closed form of the reversible drift. In the present 
section we simply quote the exact results and redirect to 
the appendices for the specific calculations. 


A. The exact reversible drift 

We have obtained in Eq. (IA24I) of Appendix the 
following exact form for the reversible drift A{x) in the 
form dm with the evidently skew-symmetric reversible 
generator 


J 


L = 


( 0 0 5^(R) 

0 0 

V -^^R) 

I- 


(43) 


The double square brackets act on arbitrary space- operation of conditional averaging and space integration, 
dependent phase functions fr{z) and denote the double 












(44) 




where ( /, 


Rpg 


is the conditional expectation ([3]) for the 


present level of description. 

The CG Hamiltonian 'H(R, p, g) is shown in Appendix 
El Eq. (lAlip to be given rigorously as 


n{R,p,g) = -ksTln 


exp|-|g^M^^ig^| \ 
(27r/det J 

.A(R,p) + $“*(R) 


Rp 


(45) 


In this expression the microscopic mass matrix is defined 
as 


N 

=^mi5^{qi) 5 jp{qi) (46) 


This matrix depends on the microscopic configuration of 
the particles and we assume that for the typical configu¬ 
rations R, p that condition the average in (l45)l are such 
that give microscopic configurations for which the inverse 
exists. 

The fluid free energy is the sum of two contributions 
T (R, p) = J'"°' (p"°') + J-“*(R, p"°i) (47) 


where the discrete solvent density is defined in Eq. 
(EU. The free energy of the solvent and the free en¬ 
ergy of interaction between nanoparticle and solvent 
are, respectively 


J^°\pso{) = 
.F“(R,psoi) = 


fcsTlnP.^'Kp.oi) 


fc^Tln 



N 

-/3^$“*(R 

i=l 



where Pgo'Kp) is the equilibrium probability that a sys¬ 
tem without the nanoparticle has a particular realization 
for the mass density. The conditional expectation 
(• • ■Y‘°^ is an equilibrium average over solvent degrees of 
freedom conditional to give the realization p^ for the dis¬ 
crete density. The fact that the free energy of the system 
in Eq. (|47ll depends on the mass density of the fluid p^ 
through the combination p™* in (1211) . which is the mass 
density of the solvent in cell p, is a non-trivial result. 


B. Approximate results for the reversible drift 

The exact but formal results (1431) , (|45p need to be ap¬ 
proximated in order to express them in terms of explicit 
functions of the relevant variables R, These re¬ 



FIG. 3: The linear for spiky approximation: The microscopic 
density field Pt[z), which is a sum of Dirac delta functions, 
each located at the particle’s position qi, is approximated 
with the linear interpolation tp^{r)p^{z) (blue line) of the 
discrete values of the density field p,j at the nodes. 


suits involve conditional expectations of the microscopic 
density fields pr('Z), gr(2)- The basic approximation that 
we will consider when computing conditional averages 
of the microscopic mass and momentum density fields is 
that these fields may be approximated by linear interpo¬ 
lations of the CG densities, this is 

Pr(z) ~ 'ip^{r)pfj,{z) 

iYz) - (49) 

A graphical representation of this approximation in ID 
is shown in Fig [31 Note that the approximation (1331) is 
equivalent to replacing the Dirac delta function (5(r — q^) 
in (fT9l) with the regularized Dirac delta function A(r, q^) 
introduced in (1331) . 

We call this approximation linear for spiky approxima¬ 
tion because Priz), as defined in Eq. m, is a sum of 
Dirac delta functions while '0^(r);O^(z) defined in (1491) is 
a piece-wise linear function of space. The approximation 
assumes that for the “typically encountered” realization 
of p, g, the above relation is well satisfied inside condi¬ 
tional expectations (• • • It is obvious that such an 

approximation makes sense only if the conditioning val¬ 
ues ppi, gfi for the densities are such that they correspond 
to a sufficiently large number of particles in cell p. Eqs. 
dMl) need to be understood in the weak sense, this is, 
valid within expressions involving space integrals. Note 
that if we multiply both sides of the approximate equa¬ 
tions (gni) with 5jy(r) and integrate over space we get an 
exact identity p^{z) = Pfj,{z) for all microscopic states 
Z', this gives us confidence in the self-consistency of this 
approximation. 

As we demonstrate in the Appendix, the linear for 
spiky approximation allows us to replace hatted func¬ 
tions with over lined functions, and to transform the dou¬ 
ble brackets |- • - into simple space averages || • • • ||. 
This transforms the exact results for the reversible drift 
into approximate but closed expressions, as we explain 
next. 
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1. Approximate mass matrix 

The microscopic mass matrix in (l46l) can be 

exactly expressed in terms of the microscopic field Pr{z) 
introduced in (fT^ . 

= \6fj,5i,p[z)\ (50) 

Note that this matrix satisfies the following exact results 

Vfj,M^u{z) = Pu[z), VuM^u[z) = Pfj,{z) (51) 

where use has been made of the first equation (1^^ . 

Under the linear for spiky approximation the 

mass matrix in (1501) becomes 

M^^{z) ~ IS^S^lilalPaiz) (52) 

and therefore, in this approximation the matrix Mp_u[z) 
depends on the microstate z only through the discrete 
density field Pa{z)- The approximation (1521) is consistent 


in the sense that it fulfills the exact properties (ISTl) . Note 
that for a function of relevant variables F{x{z)) the con¬ 
ditional expectations satisfies {F{x))^ = F{x). By using 
this property, the conditional expectation of the mass 
matrix dsni) is 

/ ^ \ Rpg _ 

~ \\Sf,6t,tj;„\\p„ = \\pSf,S4 = M^^{p) (53) 

where the interpolated mass density field p(r) is defined 
in (H^ and we have introduced the mass matrix M^^[p) 
(with dimensions of mass over volume squared) for nota- 
tional convenience. 


2. Approximate reversible generator 

In Appendix [b 1 Eq. (IB8I) . we show that under the lin¬ 
ear for spiky approximations (|49l) the exact reversible 
drift originating from the reversible operator (H51) be¬ 
comes 


/ (LR)*^'’® \ 


/ 0 

0 


\ 


/ mi \ 


( 0 ^ 



= 

0 

0 

WpsMSf^W 



an 

dpu 

-ksT 

0 

(54) 



1 -<5m(R) 





\ ^} 


1 -V“5^(R) } 



The interpolated density and velocity fields are defined 
as 

p(r) = 

g(r) = (55) 

and the double bar notation introduced in (l2^ describes 
integration over all space. The stochastic drift propor¬ 
tional to ksT emerging from the divergence of the re¬ 
versible matrix is very simple and, for the case of no 
suspended particles, indicates that the reversible dynam¬ 
ics follows a Hamiltonian dynamics, i.e., the phase space 
flow is incompressible. 


3. Approximate CG Hamiltonian 

In appendix [cl see Eq. (IC6I) . we show that under the 
linear for spiky approximation (1521) , the CG Hamiltonian 
(ggi becomes 

(R, p, g) = ^ (R, P) + (56) 


The CG Hamiltonian is the free energy of the selected 
level of description, but we refer to it as a CG Hamil¬ 
tonian because of the presence of a quadratic term in 
momenta that can be interpreted as a “kinetic energy” 
plus a “potential energy” given by the intrinsic fluid free 
energy T (R, p). This free energy is given rigorously by 

(El). 

In Appendix |Cl Eq. (IC33I) we introduce an explicit 
model for the free energy (|i7|) 

^(R,p) = 

^Peq Peq 

(57) 

where = Pfi — Peq is the density perturbation away 
from the average solvent density peq = Mt/Vt, with Vt 
being the total system volume. The motivation behind 
this model is that it gives Gaussian fluctuations for the 
solvent in the absence of any suspended nanoparticle, 
and describes in a CG manner the interaction between 
the nanoparticle and the solvent in such a way that gra¬ 
dients of density produce forces on the nanoparticle. The 
parameter cq with dimensions of speed governs the inten¬ 
sity of these forces. When the nanoparticle is simply a 
tagged solvent particle, cq = c. 
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The derivatives of the CG Hamiltonian ((551) are com¬ 
puted in Appendix [cl Eq. (|C9p 


m 

dR 

dn 

dPt, 


gjr ^^ext 

M dR 


dT 


dn 

dSp. 


= 


where the discrete velocity is defined as 


-1 




(58) 


(59) 


which is given in terms of the density dependent mass 
matrix and the momentum density field. The reason for 
introducing this somewhat involved definition for the hy¬ 
drodynamic velocity is justified by the resulting form of 
the discrete hydrodynamic equations, resembling in form 
the structure of the continuum equations. Note that in 
an “incompressible” limit in which we assume that the 
density fluctuations are very small and then = peqj 
the above expression simplifies to = Pgqgfj. because of 


^ Peq = PeqM^^. (60) 


Note that (1591) may be written as 

= = ll^^p v|| (61) 

This allows to write the interpolated momentum density 
field as 


g(r) = ^/>^(r)||5/,pv|| (62) 

If we use (1551) under an assumption of sufficiently smooth 
fields, which should apply in the limit when the grid cells 
are large and fluctuations are small, we obtain the local 
relationship 


g(r) ~ p(r)v(r) 


(63) 


proximate form for the reversible drift 


{LR)^P^ = v(R) 

{Lp,)^^^ = \\p^-V8,\\ 

= III v-V5^|| + keT^S^iR) 
dT dT 

- + 5^(R)F-‘ 

+ ^ (IIp^^V^.IIIIV'.v^II - ||p<5^Vv2||) (64) 

By conforming to the structure (El), the reversible 
drift (|64l) preserves the equilibrium distribution function 
g-pw total mass is conserved by the above 

equations, as a result of the identity (1551) . However, to¬ 
tal momentum is not exactly conserved. Since in the 
molecular ensemble ED momentum is conserved, it is 
important to conserve momentum strictly in the coarse¬ 
grained dynamics as well when F®’'* = 0, and we discuss 
this issue next. 

The rate of change of the total momentum is given by 


dPr 

dt 


dn 

dR 


llpVJ.l 


dn 

dpv 


+ ^ (llpV<5.||||^.v2|| - IpVv^l) 


(65) 


which does not necessarily vanish. The violation of mo¬ 
mentum conservation is weak, however. First, consider 
the velocity terms in (1551) . Under the assumption of 
smooth fields, Eq. (I37p applies and shows that the differ¬ 
ence of two terms in the parenthesis (last term in (1651) 1 is 
small (second order in grid spacing). Therefore, we will 
neglect the last two term in the momentum equation in 
dsi. Second, consider the terms involving the free en¬ 
ergy in (155)) . We have shown in Eqs. (IA14P and (|B9I) in 
the Appendices that the translational invariance of the 
microscopic Hamiltonian is reflected in the following ap¬ 
proximate property of the free energy 


dn ,dn „ 


( 66 ) 


which is the familiar continuum definition of velocity 
from the momentum and mass densities. In general, how¬ 
ever, (1631) does not hold identically and we prefer to define 
v(r) as the interpolant based on the discrete velocities 

d^. 


4- Approximate reversible drift 

We may perform explicitly the matrix multiplication 
in Eq. (1541) with (1551) . This leads to the following ap¬ 


relating the gradient of the free energy to the chemical 
potential . This identity implies the first two terms in 
(155)) cancel. In a way reminiscent of Noether’s theorem, 
the microscopic translation invariance (1661) implies total 
momentum conservation in Eq. (1651) . 

Unfortunately, the model for the free energy (1571) does 
not strictly respects the property (1661) . However, as we 
explain in Appendix [Cl we can restore the property (1661) 
by making the plausible approximation that the density 
field is sufficiently smooth 

Vp(R) ~ lAVpIl = J drA(R,r)Vp(r) (67) 
Recall that the reason why (EZD, which is an example 
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of (PPI) . is not an exact identity is due to the fact that 
the regularized Dirac delta is not translation invariant, 
i.e. A(r, r') ^ A(r — r'); this is the origin of the (small) 
violation of momentum conservation. If we nevertheless 
assume that the approximation (1671) is valid, then Eq. 
(IMl) is fulfilled as shown in Appendix [Cl Eq. (|B10I) and 
we restore exact momentum conservation. 

In a similar spirit, the terms involving the free energy 
in the momentum equation are computed in Appendix 
El in particular (IC38I) . with the result 

( 68 ) 

where we have introduced the “pressure” field 


standard (^ . For pair-wise interactions they are 

Lpr{z) = - V-gr(z) 

Lgriz) = -V-&y + F®’'*(qo)(5(qo - r) (71) 


where the stress tensor has the standard form 

N 


= P*v*<5(q. - r) 


2=0 


N 


9 / de5{r-(i^+ eq^) (72) 


hi=o 


Note that the stress tensor includes the nanoparticle i = 
0 in its definition. 


P(r) = 


(Cq - C^) 


2/Oe 


{p{^f - pI<^ +™ o - 


Peq 


A(R,r)p(r) 

(69) 


The time derivatives of the relevant variables 
can be obtained with (I^Ul) from the time 
derivatives of pr(z),gr(z). They are given by 


which consists of two parts, the first being the equation of 
state corresponding to the Gaussian model for the solvent 
free energy density, and the second one capturing the 
solvent-nanoparticle interaction. Note that this second 
contribution vanishes for a tagged fluid molecule, when 
Co = c. 

Inserting the result (|68p in (l64)l we get the final approx¬ 
imation of the reversible part of the momentum equation, 

= III v-+ kBTVS^(R) 
-||(5^VP||+(5^(R)F'="^ (70) 

This form exactly conserves momentum, at the expense 
of breaking the structure (HU). As a consequence, the 
equilibrium distribution that results from using the mo¬ 
mentum conserving (EOl) instead of (IM)) will be slightly 
different from oc Note that even if we has exactly 

the model of the free energy (l57)) leads to the a 
marginal equilibrium distribution of the particle position 
that is not given by the Gibbs-Boltzmann distribution 
exp {—/3$®’'‘(R)} but rather by (IG43I) . 


IV. THE IRREVERSIBLE PART OF THE 
DYNAMICS 

The dissipative matrix dS) involves the projected 
currents 6Lx = Lx{z) — {Lx)^^^\ where x{z) = 
{R, Pfi{z), g;j(z)} and LX are the time derivatives of the 
relevant variables. They are obtained by applying the Li- 
ouville operator on the position of the nanoparticle, mass 
and momentum local densities. In order to compute the 
time derivatives of the GG hydrodynamic variables it is 
useful to first consider the time derivatives of the micro¬ 
scopic local fields pr{z),gr{z) defined in (flUl) which are 


Too 

N 

2=0 

= j dvV5^{v)-(T^{z) -f F®’'*(qo)^^(qo) (73) 

The corresponding reversible part {Lx)^^^'> that is sub¬ 
tracted in the projected current has been computed in 

Eq. dMl). 

We will discuss shortly the projected current corre¬ 
sponding to the position of the colloid, which will be 
denoted by i5LR = dV. By using the linear for spiky ap¬ 
proximation (j49l) . we can approximate the time derivative 
of the density variable in (1751) as follows 

LPt^{z)- J dr tp^,{r)VS^,{r)-g^{z) (74) 

In this approximation, the time derivative of a relevant 
variable (the density) is itself given in terms of a relevant 
variable (the momentum). Therefore, the corresponding 
projected current vanishes, i.e. 5p^(z) = 0, resulting in a 
great simplification of the dissipative matrix. From Eq. 
dZSl), the projected current corresponding to the momen¬ 
tum may be expressed in the form 

5Lifj,{z)=JdrVSf,(r)-6a-r (75) 

where the fluctuations of the stress tensor are 

6a-r = d‘r(^) — (76) 

The external force term in Eq. (|73l) disappears from the 
projected current (ESI) because it is just a function of 
qo = R which is a relevant variable. 

By using (1751) . we can write the dissipative matrix T>{x) 


drV(5^(r).gr(z) 
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as a collection of Green-Kubo integrals 


1 



0 

JdrV“'^^(r')(<5V/5(0)5^“/(t) 


0 / dr'V^'S^ir') 

0 0 
\ Rpg ! a o! / \ R-Pg / r,! 

) 0 /dr/dr' (0)^6-^ W) (V“ d^(r)V^ d,(r') 


\ 


/ 

(77) 


In general, the dissipative matrix depends on the values 
of the coarse-grained variables R, p, g that condition the 
expectation values in (1771). Consider, for example, the 
colloid diffusion tensor defined as 

D(a;) = 

Indeed, even for a dilute nanocolloidal suspensions, had 
we tried to jump to the Smoluchowski level (using only 
the position of the nanocolloids as a slow variable) di¬ 
rectly, the diffusion tensor would depend strongly on the 
configuration because of the hydrodynamic interactions 
(correlations) between the particles. At our level of de¬ 
scription, however, we can assume that, to a good ap¬ 
proximation, the dissipative matrix does not depend on 
the configuration and can be approximated by its equilib¬ 
rium average, i.e., by replacing the conditional expecta¬ 
tions in (I77p with equilibrium averages. In this approxi¬ 


mation, 

DW.D-./..'P-V)D(.') (79) 

By inserting ((78l) into ((7^ and integrating over the Dirac 
delta function gives 

D“^(a;)-/ dt((5V'’(0)(5V“(t)\ (80) 

Jo ' ' 

where the average is now an ordinary equilibrium ensem¬ 
ble average rather than a constrained one. 

Under the approximation in which the dissipative ma¬ 
trix is substituted by its equilibrium average, the non¬ 
diagonal elements of the dissipative matrix dzll), which 
involve a third order tensor, will vanish because the equi¬ 
librium ensemble is isotropic and the only isotropic third 
order tensor is the null one. The dissipative matrix be¬ 
comes 




V(x) = 


0 /dr/dr'j7““'^^'v“'(5^(r)V^'(5,(r') / 
I- 


(81) 


where we have introduced a fourth order tensorial non- A. Mass diffusion 

local viscosity kernel 

( 82 ) 

The projected current corresponding to the position is 
given by 


SLR = V - = (5V 


(83) 
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where we have denoted by V = LR = ^ the velocity of 
the nanoparticle. The term is the reversible part 

of the evolution of R, given in the first equation in (IMll . 
evaluated at the microscopic value of the phase functions, 
this is 

^hydro(^) ^ 

(84) 

We expect that, being an equilibrium average, which 
is rotationally invariant, the tensor D(x) given in (1501) is, 
in fact, diagonal and of the form 

D“^ = Do5°‘f^ (85) 

Here the scalar bare dijfusion coefficient is given by 
1 r'^ 

'^^('^V(0)-(iV(t))^^ (86) 

where d is the dimensionality, and (5V is defined in (1831) 
with (|501) as the fluctuation of the velocity of the nanopar¬ 
ticle relative to the surrounding flow velocity. 

Note that the bare diffusion coefficient is different from 
the macroscopic or renormalized diffusion coefficient^ 

D = ^J\t(viO)-Vit)y'' (87) 

defined without subtracting the interpolated fluid ve¬ 
locity. We can split the renormalized diffusion coeffi¬ 
cient into two parts [a, the bare part which comes 
from under-resolved details of the dynamics occurring at 
length and time scales shorter than the ones explicitly 
represented by the discrete hydrodynamic grid, and an 
enhancement AD that comes from the advection by the 
thermal velocity fluctuations and accounts for hydrody¬ 
namic transport explicitly resolved by the discrete grid, 

D = Z?o + AD = Dq 

+ - f dt(v‘'y'^™(0)-v‘'y'^™(t) 

d Jo 

+ vhydro (0) . (t) -h <5V (0) • {t)\ (88) 

Observe that AD contains a lot of hydrodynamic in¬ 
formation because of the time lag in the time correla¬ 
tion function; during the time t hydrodynamic informa¬ 
tion (sound waves, viscous dissipation, etc.) propagates 
around the particle and affects its diffusion coefficient. 

As we elaborate in more detail in the Conclusions, the 
bare diffusion coefficient (j86|) depends on the size of the 
hydrodynamic cells, i.e., on the resolution at which hy¬ 
drodynamics is represented. By contrast, the renormal¬ 
ized diffusion coefficient (l87l) is independent of the reso¬ 
lution of the grid. However, as mentioned in the intro¬ 
duction, D is not really computable in practice in MD, 
as opposed to Do^ since the upper time limit r should be 


much larger in ((571) than in (1551) . 


B. Momentum Diffusion 

The range of the viscous kernel given in (15^ is that 
of the correlation length of the stress tensor. We will 
assume that this range is much smaller than the size of 
the cells, i.e. in the length scale in which is differ¬ 
ent from zero, the function V(5^(r) hardly changes. Note 
that the stress tensor dZH) contains the contribution of 
the colloidal particle. Therefore, a condition for this lo¬ 
cality assumption is that the colloidal particle itself is 
much smaller than the grid size. If this is the case, then 
we may adopt a local approximation 

Vrr’ !■') (89) 

and therefore the viscous contribution to the dissipative 
matrix m is 

JdrJ 

(90) 

The explicit microscopic expression for rj in (1551) is ob¬ 
tained by integrating the viscosity kernel over r, r' to get 

JdrJ ^^jyt{5crPP\0)5a--\t)y 

(91) 

where the stress tensor of the whole system is, from (1721) 

= f dr = ^pfvf + (®2) 

By using (|55|) into ([?!]) gives 

j\t(5a^^\Q)5a-'^'(t)^ (93) 

where Vt is the volume of the system. 

The viscosity tensor, being an equilibrium correlation, 
will be isotropic. The general form of the isotropic fourth 
order tensor that accounts for the symmetries of the 
stress tensor appearing in the Green-Kubo expression is 

^aa'PP' ^ ^ ^jO!/3ja'/3' ja/3'^/3a' _ 

-^^jaa'j/3/3' ('94) 

where 77 , () are shear and bulk viscosities, respectively. In 
practice, one would typically neglect the contribution of 
the nanoparticles to the viscous stress and assume that 
r]ff are the pure solvent equilibrium viscosities. 
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Finally, the dissipative matrix (|81ll becomes 


noises of the following form 


V{x)^ 


( 0 0 \ 

0 0 0 


(95) 


5af - = ^/2k^ 




+ 


ksTC 




(99) 


Note the dissipative matrix is independent of the state 
of the system due to its approximation with its equi¬ 
librium average. As a result, the stochastic drift term 
ksTOx-'Dix) in Eq. (O should be taken as zero in this 
approximation. 


where the symmetric white-noise tensor satisfies 

(Wr(i)>V()'"'(t')) = + 5 ''^' 5 ^^'''] 

X S{r — r')S{t — t') (100) 


It is straightforward to show that 


{S(Tfit)d(T'^rit')) = 2kBT5{v - v')5{t - 

( 101 ) 


C. Noise terms 


In order to construct the Ito SDE @ for the present 
level of description, we need to specify the noise terms 
The variance of the noise is given by the 
Fluctuation-Dissipation balance where the matrix 
D(x) is given by ([Ml) . From the structure of this matrix 
we may infer that = 0 and 

(^W^(o) = 2fcBrDo<5(t-0 


and, therefore, the correlation of the random stress is a 
white noise in space and time, proportional to the viscos¬ 
ity tensor. Now we use the following expression for the 
piece-wise constant gradient of the finite element linear 
basis functions (Hi 

VV'.(r) =^be,0e.(r) (102) 

Gi/ 

where labels each of the sub-elements of the node v, 
is a constant vector within the sub-element that 

is pointing towards the node ly and 6f.^ (r) is the charac¬ 
teristic function of the sub-element e^. 


^^(O^(i')) = 2fcBrj7““''5^'||V“'(5^V^'(5,||5(f - t') 

(96) 

We need to produce next explicit linear combinations of 
white noise that give rise to the above variances. While 
the velocity noise term is very simple 

(7R 

— {t) = ^2kBTDoWit) (97) 

at 


The projected current, can be written, therefore, as 

[ dr'^h^_^9e,{r)dcr‘^^ (t) (103) 

By using the model (IMl) for the projected stress tensor 
and equating the random term with the projected 
current SLg^ we have the following explicit model for 
the random forces 


where W(t) is a white noise, the explicit form of the 
random force is not so obvious and will be considered 

at 

next. 


= (104) 

et/ 


The noise term in the theory of CG is just a modelling 
of the projected current appearing in the Green-Kubo 
expression as a white noise. For this reason, it is 
useful to look at the structure of the projected current in 
Eq. (HSl) 


where the random stress tensor of the sub-element e^, is 
given by 


(t) = y/2k^ 


(<) - EW 


I (r)5^y (98) 





(105) 


We will model 


as a linear combination of white 


Here, we have introduced a symmetric matrix of white 
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noise processes associated to each sub-element 

wr(<) = J drde(r)wr(0 (106) 

These symmetric white-noise processes are independent 
among elements due to (llOOl) 

= See'[s^^'+ s''f^'s>^''']s{t -1') 

(107) 

The noise term (11041) is a discrete divergence of a dis¬ 
crete random stress tensor. The discrete random stress 
tensor is an independent stochastic process associ¬ 
ated to each sub-element. It is a matter of calculation 
to check that the postulated noise term in (I104I105P 
with the white noise per element We^(t) satisfying (I107p . 
gives precisely the FD balance in (IM)) . 

Note that the noise (I104p contains the matrix 
which is the inverse of defined in (I^Sl) . The elements 
ofM^, are proportional to the volume (area in 2D) of the 
overlaping region between two hydrodynamic cells which, 
in turn, scales as the typical volume of the hydrodynamic 
cells. Therefore, the stochastic force -^{t) scales with 
the inverse of the square root of the cell size. Larger cells 
are subject to smaller fluctuations, in accordance with 
the usual concepts in equilibrium statistical mechanics. 


V. FINAL APPROXIMATE DYNAMIC 
EQUATIONS 

We now have all the ingredients to construct the SDE 
(|9|) for the chosen coarse-grained level of description. By 
collecting the reversible part (l64)) with (|70)) and irre¬ 
versible part of the dynamics given by 'D-dx'H (where 
the dissipative matrix is (1951) and the derivatives of the 
CG Hamiltonian are in ((58l) l. the final SODEs for the 
selected CG variables are 

^ = WRi - ^ 

dt ^ ’ keT OR keT dt 

^ = ||g V-Vd^l + fcBrV<5^(R) - Ild^VPlI + d^(R)F‘^"' 

+ + (i + 0 ^ 

These equations are the main result of this paper. Recall 
that the double bar denotes the spatial average defined 
in (|2^ . and the overlined symbols denote interpolated 
fields out of the discrete values as in, for example, v(r) = 
V[^'0[^(r), etc. The velocity is given in terms of 
in (15^ . The pressure equation of state is given in (1551) 


and the gradient of the free energy (1571) is given by 

(109) 

oR Peq Peq 

see (1571) for the definition of the notation ||AVp||. The 
SDEs (11081) are closed and explicit in the relevant vari¬ 
ables. 


A. Physical meaning of the different terms in the 
dynamic equations 

The first equation in (11081) governs the evolution of 
the position of the nanoparticle. The first term v(R) is 
purely reversible and says that the nanoparticle follows 
the interpolated velocity field of the fluid. This is a purely 
kinematic effect due to the fact that the momentum of 
the fluid contains the contribution due to the nanopar¬ 
ticle. It has nothing to do with any force that the fluid 
may perform on the particle which are described by the 
second contribution. This contribution is proportional to 
the bare mobility Dq/ksT, given in terms of the bare dif¬ 
fusion coefficient Dq introduced in (1551) through a Green- 
Kubo relation. This term involves the (minus) gradient 
of the free energy P(R, p), which plays the role of a po¬ 
tential of mean force for the nanoparticle given explicitly 
in (I109p . As seen in (IG30p the force due to the fluid 
on the nanoparticle involves the gradients of the solvent 
density. The presence of the two parameters cq, that 
is due entirely to interactions of the nanoparticle with 
the solvent particles, and c, which is due to interactions 
of the solvent particles with themselves alone, indicates 
that — is not simply the force that the solvent ex¬ 
erts on the nanoparticle. Note that in the limit when 
the nanoparticle becomes just a tagged solvent particle, 
which is realized for cq —>■ c, — given in (11091) vanishes. 
The third term in the position equation in (llOSp is due 
to the external force that obviously affects the motion of 
the nanoparticle. Finally, the nanoparticle is subject to 
an explicit noise term ^ whose variance is given by the 
fluctuation-dissipation balance relation (1951) . This term 
will produce Brownian motion of the nanoparticle, in ad¬ 
dition to the advection by the fluctuating velocity field 
v(R). In order to not “double count” the noise in the 
Brownian motion of the particle m , the diffusion coef¬ 
ficient that governs the amplitude of the random noise 
term ^ is given in terms of the bare diffusion coeffi¬ 
cient Dq in (I86p . and not by the renormalized diffusion 
coefficient D defined in (l87P . 

The second equation in (|108l) gives the evolution for the 
discrete mass density and has the form of a discrete 
continuity equation. This evolution is purely reversible 
due to the fact that, very approximately, the time deriva¬ 
tive of the mass density is given in terms of the momen¬ 
tum density, which is a relevant variable. Therefore, the 
projected current vanishes and so do the Green-Kubo 
coefficients, i.e., there are no Brenner diffusion terms, as 
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argued in (5^. 

The third equation in (jl08l) governs the discrete mo¬ 
mentum density g^. It has the structure of a discrete 
version of the fluctuating isothermal compressible Navier- 
Stokes equations with some modifications due to the 
interactions with the nanoparticle. The first term in 
the momentum equation is a convective non-linear term 
quadratic in the discrete momenta, which corresponds 
to the usual convective term in the Navier-Stokes equa¬ 
tions. The second term originates from the stochastic 
drift ksTdx-L term and can be interpreted as an osmotic 
pressure term due to the presence of the nanoparticle. 
The third term is reminiscent of the pressure gradient 
term in the usual Navier-Stokes equations. The pres¬ 
sure equation of state is given by the pressure due to the 
Gaussian model for the solvent, plus a pressure correc¬ 
tion term (proportional to the difference of the squares 
of the speeds of sound) that describes the interaction 
between the solvent and the nanoparticle. Finally, the 
term proportional to F®’'* in (11081) describes the effect 
that, because the discrete momentum variable contains 
the contribution due to the nanoparticle, any external 
force on the nanoparticle will translate into a force on the 
fluid itself. All the terms discussed so far in the momen¬ 
tum equation are purely reversible. The only irreversible 
terms in the momentum equation are proportional to 
the viscosities rj, and correspond to the usual viscous 
terms involving second space derivatives in the Navier- 
Stokes equations. Finally, the term dg^ is the random 
forces with explicit form given in (11051) and whose ampli¬ 
tudes are dictated by the fluctuation-dissipation balance 
in (IMl)- 

Note that when cq = c a number of terms in the equa¬ 
tions above drop out and the equations simplify consider¬ 
ably. This happens, for example, when the distinguished 
particle is simply a tagged fluid molecule. This may also 
be a good approximation for neutrally buoyant particles 
that do not have a strong chemical interaction with the 
surrounding fluid, and the majority of prior work in the 
literature has in fact used the simplified model cq = c, 
with the notable exception of (3^ . 


B. Scope and general properties of the dynamic 
equations 

The validity of the SODEs (I108D is limited to situa¬ 
tions in which the values , g^ of the relevant variables 
are such that give a large number of solvent particles per 
hydrodynamic cell and, at the same time, give values that 
do not differ very much from one cell to its neighbors. In 
other words, the interpolated fields p(r),g(r) need to be 
smooth on the hydrodynamic cell length scale. These 
assumptions imply that the validity of the equations is 
restricted to situations in which thermal fluctuations are 
small. Correspondingly, we have assumed that the sol¬ 
vent density fluctuations are Gaussian. This precludes 
the study of other interesting phenomenology like liquid- 


vapor phase transitions, for example. However, it is a 
sufficiently simple and physically realistic model in many 
situations of interest. Concerning the nanoparticle, it is 
assumed that it is smaller than the hydrodynamic cell 
and it is, therefore, a subgrid nanoparticle. 

The SODE (jl08l) conserve exactly the total mass of the 
system defined in (1551) . In the absence of external forces 
acting on the nanoparticle, = 0, the total momen¬ 
tum is also exactly conserved by the equations. This is 
just a reflection of the definition (ITTI) of the discrete mass 
and momentum “fields” in terms of the basis functions 
that satisfy the partition of unity property (1551) . Momen¬ 
tum conservation is a direct consequence of translational 
invariance and we restored exact momentum conserva¬ 
tion in our approximate equations by restoring transla¬ 
tional invariance of our free-energy model. 

Discrete fluctuation-dissipation balance (DFDB) is a 
crucial property that has been carefully maintained in 
prior work that relied on phenomenological equations, 
see for example [s^ or Appendix B of [Sg. A key com¬ 
ponent of DFDB is the energy conservation property that 
any work done by the external forces on the suspended 
nanoparticle must be converted exactly into kinetic en¬ 
ergy of the fluid. In the terminology of Refs. [2^ [s^ [s^ , 
this means that the linear operator (matrix) used to in¬ 
terpolate the (discrete) fluid velocity to the particle, as 
represented by the term (i?) in the first equation of 
(11081) . is the adjoint of the linear operator used to spread 
the force applied on the particle to the fluid, as repre¬ 
sented by the term (R) (R) in the last equation 

in (11081) . This energy conservation follows directly from 
the skew-symmetry of the reversible operator (1541) . 

If the reversible drift were exactly in the form Ld^H — 
ksTOx-L^ it would automatically maintain DFDB, this is, 
the equilibrium distribution function would be oc e~^^. 
However, the smoothness approximation taken in order 
to arrive at the model (jl08p imply that this is true up 
to small second order terms in the lattice spacing. Our 
selection of the model for the free energy is not exactly 
translational invariant, i.e. it does not satisfy (1661) ex¬ 
actly. If it were, as shown in Appendix [Cl then we would 
obtain that the marginal distribution P®‘i(R) for the po¬ 
sition of the particle would be given exactly by the baro¬ 
metric law (reduced Gibbs-Boltzmann distribution), 

P®®'(R) - exp {-/3$®’^fyR)} (110) 

However, the violation of translation invariance implies 
that the resulting probability distribution is given by 
(|G43p instead, and the true barometric distribution is 
obtained only in the incompressible limit c —>■ oo or if 
Co = c (e.g., a tagged particle). 
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C. The continuum equations 


VI. DISCUSSION AND CONCLUSIONS 


We have obtained the SODEs (11081) from the Theory 
of Coarse-Graining. It can be shown that the same equa¬ 
tions can be obtained from a Petrov-Galerkin discretiza¬ 
tion (see [131 for an illustration using the same basis func¬ 
tions as used here) of the following system of stochastic 
partial differential equations (SPDEs) 

= J (irA(r, R)v(r) 

J- Peq J 

, Dp pext , ^ 

keT dt 

dtp{r,t) = -V-g 

ftg(r,t) = -V-(gv) - fcsrVA(r,R) 

- VP(r) -f r®’'*A(r,R) 

TryV^v-l-(|-fC) V(V-v) +(111) 
where v = g/p, and the pressure is given by 

^ {p{rf - Peq) + —^A(R,r)p(r) 

2peq ^ Peq 

( 112 ) 

The random velocity dH/dt is given in (l97l) . and the ran¬ 
dom stress tensor is given in (IMl) . The equations 
(11111) are very closely related tophenomenological equa¬ 
tions used in prior work with some differences 

that we further discuss in the Gonclusions. 

The Petrov-Galerkin method in its most pedestrian 
form has three steps: 1) Multiply the equations (11111) for 
the hydrodynamic fields with the basis functions (5^(r) 
and integrate with respect to space. 2) Define the dis¬ 
crete variables P/j = / dr(5^(r)p(r), etc. 3) Approximate 
the fields in the right hand side of the equations (11111) 
with the linear interpolations p(r) = '0^(r)p^, etc. This 
procedure applied to (11111) then leads to (11081) . As an 
example, let us consider the hrst term in the equation 
of motion for the particle, representing the advection by 
the fluid velocity. Replacing the velocity with its linear 
interpolant we get 

y V (r) A (r, R) dr^ J (r) A (r, R) dr 

= J '(IJr. (r) A (r, R) dr 

= W = V (R), (113) 

where we used the property (1351) . The right hand side is 
exactly our discretization (derived here from the micro¬ 
scopic dynamics!) of the term on the left hand side. The 
rest of the terms are discretized in a similar way. 


By performing a systematic coarse-graining procedure 
based on the Zwanzing projection operator, we have de¬ 
rived a system of stochastic ordinary differential equa¬ 
tions (I108D describing the dynamics of a nano-sized par¬ 
ticle immersed in a simple liquid. A key to the proce¬ 
dure was the use of a dual set of linear basis functions 
familiar from finite element methods (FEM) as a way to 
coarse-grain the microscopic degrees of freedom. Another 
key ingredient was the use of a “linear for spiky” weak 
approximation which replaces microscopic “fields”, i.e., 
sums of delta functions centered at the fluid molecules, 
with a linear interpolant in the FEM basis set. These 
two steps enabled us to obtain closed approximations for 
all of the terms in the reversible or non-dissipative dy¬ 
namics, in a manner which gives them a clear physical 
interpretation and preserves the correct structure of the 
equations. Notably, the reversible dynamics preserves a 
discrete Gibbs-Boltzmann distribution to high accuracy. 
For the irreversible or dissipative dynamics, we approx¬ 
imated the constrained Green-Kubo expressions for the 
dissipation coefficients with their equilibrium averages, 
and assumed a local form for the viscous dissipation suit¬ 
able when the hydrodynamic cells contain a large number 
of fluid molecules. 

The coarse-grained equations we derived here can be 
seen as a particular Petrov-Galerkin FEM discretization 
of a system of continuum stochastic partial differential 
equations (llll|l coupling the familiar isothermal fluctu¬ 
ating Navier-Stokes (FNS) equations with the Brownian 
motion of the immersed particle. These equations are 
similar in structure to phenomenolo gical eq uations used 
in a number of prior works [TlL [H, 1^ . l34l . , and 

therefore provide a justification for those types of mod¬ 
els via the Theory of Goarse Graining. This is not just 
an academic exercise, but one that also has some im¬ 
portant practical utility. First, our derivation provides 
Green-Kubo expressions for transport coefficients, no¬ 
tably, for the hare diffusion coefficient which was phe¬ 
nomenologically added in [s^ as a way to account for 
under-resolved microscopic details that cannot be cap¬ 
tured with a hydrodynamic approach. Our derivation 
also introduces novel terms that come from the micro¬ 
scopic interaction between the suspended particle and the 
liquid molecules; these terms give additional modeling ca¬ 
pability to account for more microscopic information in 
the coarse-grained description. Another, perhaps unex¬ 
pected, benefit of the microscopic derivation was that it 
lead directly to a discrete form of the divergence of the 
stochastic stress which obeys the fluctuation-dissipation 
balance relation. In more empirical approaches, such a 
structure has to be either guessed, or constructed from 
suitable discrete stochastic fluxes and a pair of discrete 
divergence and gradient operators that are skew adjoints 
of one another (60|. 

The coarse-graining procedure carried out here can be 
viewed as a systematic derivation of the isothermal FNS 
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equations from molecular dynamics. Formal derivations 
of these equations have been done before many times, 
see for example early work including non-linearities in 
[13, mi , however, these derivations lead either to lin¬ 
earized equations or to ill-defined nonlinear SPDEs ex¬ 
hibiting an ultraviolet catastrophe. As we argued in more 
detail in (^ . the proper way to interpret such formal 
nonlinear SPDEs is to first discretize them by applying a 
systematic discretization procedure, for example, using a 
Petrov-Galerkin weak formulation. The justification for 
this prescription is the fact that here we obtain exactly 
the same set of discretized SPDEs by systematic coarse 
graining. This gives a direct link between the “bottom- 
up” approach of going from microscopic to mesoscopic 
equations, and the “top-down” approach in which one 
starts from continuum PDEs and formally adds white- 
noise stochastic forcing and then applies a standard com¬ 
putational fluid dynamics (CFD) method to the resulting 
equations. Renormalization techniques should be applied 
to the discrete equations rather than the continuum ones 
in order to systematically increase the coarse-graining 
scale from the mesoscopic to the macroscopic in order 
to recover a continuum limit, where predictions of physi¬ 
cal quantities like space-time correlations do not depend 
on the lattice spacing h as h —0. 

The crucial link between the top-down and bottom-up 
approaches was already foreseen in [i^, and then ex¬ 
plicitly demonstrated on a significantly simpler micro¬ 
scopic model in the Ph.D. thesis [s^. At the same time, 
the new derivation given here significantly improves on 
the earlier derivation which did not consider a sus¬ 
pended nanoparticle, in three key ways. Firstly, here we 
account for the presence of a suspended particle. Sec¬ 
ondly, by using a dual set of basis functions, the result¬ 
ing discretization is second-order rather than first-order 
accurate as the earlier derivation based on a single set 
of basis functions . Thirdly, in the present work the 
discrete equations (I108D have a very precise relation to 
the continuum equations (11111) . rather than simply being 
reminiscent of some “sensible” discretization. We have 
given an explicit prescription of how to connect the two 
worlds of MD and CFD: use the same dual set of basis 
functions when coarse-graining as you do when discretiz¬ 
ing. We believe this prescription can be applied to a 
variety of other problems, however, as this work shows, 
the bottom-up approach requires a lot more work to com¬ 
plete than the top-down approach. 


A. Relation to phenomenological models 

The continuum equations (11111) we proposed here bear 
a strong similarity, but also some crucial differences, with 
existing models. In order to explain the relation to prior 
work more clearly, let us review a variety of existing mod¬ 
els starting from more “refined” to more “coarse.” A more 
formal mathematical presentation of these levels of de¬ 
scription is given by Atzberger , here we give a phys¬ 


ical summary. In many works an incompressible approx¬ 
imations is made in order to eliminate fast sound waves 
from the model [s^, [s^ ■ If one is interested only in the 
long-time diffusive dynamics the fluid inertia can also be 
eliminated by taking an overdamped limit [Til. 1^1^ . 
Here we focus on the coupling between the nanocolloidal 
particle and the fluctuating fluid and not on the specifics 
of the fluid equations. 

In a number of “point particle” frictional coupling ap¬ 
proaches, as reviewed in detail in [dlj, the colloid’s ve¬ 
locity is also included as a physical variable and a phe¬ 
nomenological “friction” force proportional to the velocity 
of the colloid relative to the local fluid velocity is added 
to an inertial particle equation. As we discussed in more 
detail in Section IIIBi such a level of description is not 
suitable under our assumption that the colloid is smaller 
than the typical size of the hydrodynamic cells. Instead, 
we include the momentum of colloid in the total hydro- 
dynamic momentum field. We note that, although not 
usually presented in this way, the same is actually true to 
a large extent for the model used in [^ because what is 
called the mass of the colloid is actually the excess mass 
of the colloid over the expelled fluid [sj, IHl . This is be¬ 
cause the inertia of the expelled fluid is already included 
in the FNS equations, which are assumed to apply ev¬ 
erywhere including the volume occupied by the particle. 
Therefore, part of the momentum of the particle is in fact 
included in the “fluid momentum”. For this reason, we 
find it difficult to imagine how one can justify the fric¬ 
tional point particle coupling model from a microscopic 
derivation. 

In 1^ . [s^ l , an instantaneous inertial coupling is pro¬ 
posed in which the particle is forced to follow the local 
fluid velocity; the main difference is that [s^l considers a 
compressible fluid, while l38ll focuses on the incompress¬ 
ible limit. As shown in |62l |. the instantaneous coupling 
IHisi can be derived as a limit of the frictional cou¬ 
pling when the friction coefficient becomes very large. In 
the language of the TCG, the (fast) momentum of the 
particle is no longer included as a relevant variable, in¬ 
stead, just as in our description, a total momentum field 
is defined [s^ . This total momentum field follows an 

equation which has a similar structure to the momentum 
equation in (Hill) , see for example Eq. (16) in [Hi in 
the incompressible limit. It is important to note that, 
even in the limit of an incompressible liquid, the density 
p appearing in Hill) is not constant, rather, it includes 
the contribution from the colloid. Therefore, for a dense 
particle (e.g., gold nanocolloid) the discrete density will 
be larger at the nodes in the vicinity of the particle. This 
“excess inertia” is explicitly included in the model in Refs. 
[H.[H| (see for example the left-hand side of (13) in S) 
by “spreading” the excess inertia to the fluid grid. In 
the equations derived here the excess inertia is hidden 
in the definition of the coarse-grained density to include 
the contribution from the colloid, and the fact that the 
equation of state is only applied to the solvent part of 
the density 
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The terms in the particle equation related to the bare 
diffusion coefficient Dq do not appear in either or 
M,[62|; these terms are suggested in Appendix B.l of 
38l | but not included in numerical simulations. A bare 
diffusion coefficient is also introduced in the theoretical 
work [llj but it is argued there that this term should 
somehow be small. Interestingly, a renormalization of 
the diffusion coefficient similar in spirit to Dq is present 
in the frictional coupling formulation and can be ex¬ 
pressed in the Einstein form Dq = ksT /y, where 7 is the 
phenomenological friction coefficient, see Eq. (290) in 
|4l|. In the limit of infinite friction, which is how (s^ . 1^ 
derives the instantaneous coupling equations, Dq —>■ 0. 
However, this is “throwing the baby out with the bath¬ 
water” and is not consistent with our microscopic deriva¬ 
tion. The fact that Dq > 0 is easy to appreciate: the 
sum D = Dq -|- ad is a physical parameter that can be 
measured and is independent of the grid spacing (i.e., 
the coarse-graining length scale) while AD < D depends 
strongly on the grid spacing, as we explain in more de¬ 
tail shortly. In our derivation, Dq emerges naturally as 
does a Green-Kubo expression for it, giving it a precise 
microscopic interpretation. The fact our instantaneous 
coupling equations with bare diffusion cannot be consis¬ 
tently derived from the frictional coupling formulation 
[dll points to the lack of a microscopic foundation of 
that formulation, and justifies once again the advantage 
of systematic bottom-up approaches over phenomenolog¬ 
ical ones. 

The stochastic thermal drift term —fenTyAfr, R) in 
the momentum equation is (wrongly) missing in |34| ; the 
term is also (rightfully) missing in the frictional coupling 
[4l| formulation. This term ought be there for instan¬ 
taneous coupling, as explained in Appendix B of 
based on fluctuation-dissipation balance arguments, and 
derived by taking the limit of infinite friction in the fric¬ 
tional coupling in (^ . This osmotic pressure contribu¬ 
tion from the particle, spread to the fluid via the reg¬ 
ularized delta function, can be seen as coming from the 
eliminated fluctuations of the particle velocity around the 
local fluid velocity [s^ [^. In the inertial coupling for¬ 
mulation, as explained in Appendix B of and also 
in 0 , this osmotic pressure is split into two pieces just 
like the particle momentum is split into two pieces, one 
piece attached to the fluid momentum and another excess 
piece. When the two pieces are added together one cor¬ 
rectly recovers exactly the osmotic pressure term given 

incmna. 

The gradient of the solvent pressure appears in all 
phenomenological models, and seems very natural but it 
should be recognized that this is only an approximation. 
Notably, the approximation consists in the assumption 
that this pressure is given by the equation of state of the 
fluid in the absence of the colloidal particle. This ap¬ 
proximation is exact for a labeled or tagged particle of 
the fluid, i.e., in case of self-diffusion. One can argue that 
the same should approximately hold for colloidal particles 
that have a similar structure to the fluid, notably, that 


have the same density and compressibility as the fluid. 
Balboa et al. have proposed adding an additional 
excess pressure term to account for a different compress¬ 
ibility of the colloid relative to the surrounding fluid, see 
(I116p . These terms are postulated on a phenomenological 
basis. 

In this work we derived equations containing similar 
terms, however, as already explained, our model for the 
free energy differs from that used in and our final 
equations (Hill) are different from their, which can be 
written in our notation as 


= 

dt 

J A{R — r)v(r) dr, 

(114) 

dtp = 

-V-g, 

(115) 

dtg = 

-V-{gv)-{kBT)VA{R-r) 



-VP + r®’^*A(JR -r) + v-cr. 

(116) 


where 

P = c2(p(r)-peq) 

+ V{cl - c^)A{R -v) f A{R - r')(p(r') - Peq) dr' 

(117) 

where the regularized delta function A{R, r) = A(il — r) 
is given by a Gaussian-like isotropic kernel A(r) of width 
comparable to the hydrodynamic radius of the nanopar¬ 
ticle and that integrates to unity. Here the fluid velocity 
is defined via v = g/p, and a denotes the viscous stress 
(deterministic and fluctuating components) with a form 
identical to ours. The volume associated to the particle 
can be expressed in our notation as V = mo/Peq but is 
interpreted differently in (s^l to be a geometric rather 
than an inertial quantity. Note that we have set here the 
excess mass of the particle [s^l over the fluid mg = 0 
because in our notation p includes the total mass of the 
particle (see additional discussion in Section IVTl) . 

The similarities between our formulation and the for¬ 
mulation of Balboa et al. are evident, but there are also 
some notable differences. The equations HI6I) do not in¬ 
clude bare diffusion (i.e., Dq = 0) and therefore a number 
of terms are missing from the particle equation. In the 
model of Ref. [s^l the parameter cq is used to represent 
the speed of sound (i.e., the isothermal compressibility) 
inside the colloidal particle, while in our model cq models 
the mean force that the colloid experiences in a density 
gradient, see (IC30I) . Balboa et al. justify the “parti¬ 
cle compressibility” pressure contribution proportional to 
Cq — starting from a quadratic contribution to the free 
energy density of the form 

•^comp = ^ (/ ~ ~ 

whereas the “continuum” analog of our linear model for 
the interaction free energy found in the last term in (1571) 
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is 


Tint = V{cl 


j A{R - r')(p(r') - po) dr' 


D Ri Dse 


keT 

arjR’ 


It is important to emphasize, however, that our discus¬ 
sion of similarities to prior work above only concerns the 
formal continuum formulation (jllll) and focuses on the 
structure of the equations and the physics of the various 
terms. Our fully discrete formulation (|108l) is completely 
new and is different in many crucial ways from existing 
discretizations. The first difference is that the discretiza¬ 
tion of the FNS equations is based on a second-order con¬ 
servative FEM method, rather than the more comm only 
used finite-difference or finite-volume approach [^ . 
A second difference is that in all prior models we are 
aware of, the regularized delta function is used to repre¬ 
sent the particle itself, and its width is chosen to be on 
the order of the hydrodynamic radius of the particle. As 
such, the regularized delta function is attached to the par¬ 
ticle, i.e., it is taken to be a smeared delta function kernel 
such as a Gaussian kernel or the tensor-product kernel of 
Peskin centered at the particle position. By contrast, 
in the formulation (Iflfl) the regularized delta function 
represents the coarsening of the solvent dynamics and its 
width is assumed to be larger than the subgrid colloidal 
particle. As such, our discrete delta function kernel is 
attached to the grid and is not centered around the par¬ 
ticle; this is most obvious in (I108D where it is clearly seen 
that in the equation for the (fixed in space) grid node p 
the regularization enters via the basis function (5^ associ¬ 
ated to that node. Note that in the immersed-boundary 
formulation of Peskin used in a number of prior works 
[ 41 II the width of the discrete delta function 
is tied to the grid spacing just as it is for our regularized 
delta function, however, in those prior works the regular¬ 
ized kernel is still centered around the particle. This was 
proposed by Peskin as a very effective way to maximize 
translational invariance of the particle-grid interactions 
we expect our formulation will not perform as well 
in terms of translational invariance because it was not 
explicitly constructed with that goal in mind. 


B. Renormalization of the diffusion coefficient 

Let us consider, for a moment, a single freely-diffusing 
isolated spherical nanoparticle suspended in a quiescent 
fluid. At large time scales, the particle will perform a 
standard Brownian motion with a renormalized diffusion 
coefficient D given in (l87ll by the familiar Green-Kubo 
integral of the particle’s velocity autocorrelation func¬ 
tion. Since this quantity only involves the nanoparticle 
position, it cannot depend on how we chose to perform 
the coarse graining of the fluid. In particular, it must 
be a number independent of the typical grid spacing h. 
For sufficiently large Schmidt numbers |33| we expect it 
to be well-predicted by the Stokes-Einstein formula (in 


where R is the radius of the spherical particle and a is 
a coefficient that depends on the boundary conditions 
applicable at the surface of the sphere, equal to Gtt for 
a stick surface and 47r for a slip surface, or something 
in-between for more realistic models [M , [63 • 

In the fluctuating hydrodynamic formalism presented 
here, the renormalized diffusion coefficient D is split into 
a bare part Dg and a renormalization AD defined by (l88)) . 
Let us try to get a more quantitative understanding of the 
diffusion enhancement AD for our specific discretization 
(approximation) of the equations. First, we can replace 
the instantaneous interpolated fluid velocity by 

its approximation v(R) = v^^^(R). Second, we can 
ignore the cross-correlation terms since V evolves on a 
much faster time scale than the hydrodynamic fields and 
can be assumed to be a white-noise process uncorrelated 
with v(R). This gives the approximation to the diffusion 
enhancement produced by our discrete equations, 

AD{R) = dt{v{Rm-^R{t))r 

(118) 

In the overdamped limit of large Schmidt numbers (see 
[ 3 ^ for corrections at moderate Schmidt numbers), the 
particle moves much slower than the hydrodynamic cor¬ 
relations decay, and one can express the diffusion en¬ 
hancement in terms of the equilibrium correlation of the 
fluid velocity (conditional on the particle being fixed at 
a particular location), 

ADiR) = -j dt ^^(R) {v^iO)-v^,it)r^tP^,iR), 

(119) 

which can in principle be computed exactly by linearizing 
the fluid equations around a quiescent state. Note that 
AD depends on R explicitly; for confined systems this 
dependence is physical but for translationally invariant 
system such dependence is a discretization artifact that 
is hopefully small. Note that the immersed-boundary 
discrete delta function used in Refs. [ 2 ^, [H, is 

specifically designed to obtain such translational invari¬ 
ance on regular grids to a high accuracy (^ . 

We can obtain a physical estimate for the diffusion due 
to advection by the thermal velocity fluctuations by as¬ 
suming that the discrete velocities are consistent with a 
Petrov-Galerkin procedure applied to continuum equa¬ 
tions. This allows us to approximate the discrete veloc¬ 
ity v^ in terms of a continuum fluctuating field v(r, t) as 
v^(t) = f drSf^(r)v(r,t). If we substitute this in (III9I) 
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we obtain the estimate 


ADCR) = 


X 


X 


J drdr' ^/>^(R)(5^(r) 

J drdr' A(r, R) 

Q^°“dt(v(r,0)-v(r',t)r) A(r',R). 

( 120 ) 


If one assumes that the evolution of v(r, t) can be de¬ 
scribed by a fluctuating Stokes equation (i.e., linearized 
incompressible flow), the time integral can easily be ex¬ 
pressed in terms of the inverse Stokes operator (i.e., the 
Green’s function for Stokes flow) [Hi. In this case the 
relation (I120p can directly be matched with Eq. (10) in 
[Hi (see also Eq. (288) in 111), where the regularized 
delta function is denoted with A(r, R) —>• (5a(r — R). In 
our notation Eq. (10) in [Hi becomes, 

AZ)(R) = —— [ drdr' A(r,R) 

V J 

)A(r',R), (121) 


-Trace G (r, r') 


where G is the Green’s function for Stokes flow (Oseen 
tensor for an unbounded domain at rest at infinity). This 
is nothing else but an Einstein formula relating the dif¬ 
fusion coefficient with the mobility of the particle, i.e., 
with the linear response of the particle to a weak applied 
force. 

A simple calculation based on the expression for G in 
Fourier space, or, equivalently, based on replacing the 
continuum Green’s function with its discrete equivalent, 
estimates that in three dimensions [Hi 


AD 


keT 
a'rjh ’ 


where h is the width of the regularized Delta function 
(i.e., the grid spacing), and a' is a coefficient that de¬ 
pends on the geometric details of the grid. This suggests 
that 


which must be non-negative, i.e., it must be that aR < 
a'h, which is consistent with the assumption that the 
nanoparticle is smaller than a typical grid cell. Observe 
that at large Schmidt numbers AD can be expressed 
purely in terms of geometric quantities and the equi¬ 
librium (discrete) fluid correlation functions, and should 
therefore depend mildly if at all on the details of the 
interaction between the particle and the fluid. This sug¬ 
gests that it is the bare diffusion coefficient that must 


capture essentially all of the microscopic details such as 
slip versus no-slip on the particle surface or layering of 
the fluid around the particle. 

In the microscopic derivation presented here, the bare 
diffusion coefficient Dq is to be computed using (I5S1) from 
the Green-Kubo integral of the autocorrelation function 
of the particle peculiar velocity, i.e., the velocity relative 
to the local (interpolated) fluid velocity. Only a combi¬ 
nation of molecular dynamics and fluctuating FEM cal¬ 
culations can tell us whether the effective diffusion coeffi¬ 
cient D = Dq + AD indeed be a constant (approximately) 
independent of the grid resolution. In prior work based 
on p henomenological fluctuating hydrodynamics theories 
[n|. Do was treated as an adjustable parameter that is 
chosen so as to give a desired (input) effective D, since 
AD follows from the discretization of the fluid equations 
and cannot be adjusted independently. This is similar to 
how one can treat the fluid-particle interaction strength 
parameter cq as fitting parameters used to match the 
coarse-grained and particle dynamics as best as possible, 
instead of computing it from its (approximate!) micro¬ 
scopic definition (IG23jl . 


C. Future Directions 

The work described here is purely theoretical and pro¬ 
poses a model with the correct structure but leaves a 
number of terms to be approximated and modeled. As 
such, the usefulness and accuracy of our equations can¬ 
not be judged until a number of numerical studies are 
performed. 

Firstly, a temporal discretization needs to be devel¬ 
oped to go with the spatial discretization (I108D : the re¬ 
quired tools are readily available mill- Secondly, it is 
important to study the numerical aspects of the Petrov- 
Galerkin FEM discretization developed here using exist¬ 
ing numerical analysis tools (H . fH | and compare to exist¬ 
ing discretizations. Lastly, one needs to include immersed 
particles in the numerics as well and study a number of 
standard test problems to evaluate the performance of 
(I108|) as a standalone method for simulating dilute col¬ 
loidal suspensions. Particular emphasis should be payed 
to the violations of discrete fluctuation dissipation bal¬ 
ance and of translational invariance, both of which are in 
a formal sense second order in the grid spacing, but may 
be significant in practice at scales comparable to the grid 
spacing. 

It is an important task for future work to perform 
molecular dynamics (MD) simulations and compare the 
results to the coarse-grained description proposed here. 
We expect that if the grid cells are too small we will see 
unphysical artifacts, and if the grid cells are too large, the 
MD simulations will become unfeasible. By confirming 
whether the correct effective (renormalized) diffusion co¬ 
efficient is obtained over a reasonable range of grid spac- 
ings, we can access how good the approximations made in 
our coarse graining theory are, and ultimately how useful 
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the proposed equations are in practice. One should begin 
these studies with a single nanocolloidal particle in sol¬ 
vent in a periodic box of varying sizes, perhaps starting 
in two dimensions, where there are very strong (in fact, 
asymptotically dominant rather than decaying) finite size 
effects [TH - 

A key question that we did not fully address here is 
how to compute the various coefhcients that appear in 
(11081) . We already discussed the subtlety of this issue 
for the bare diffusion coefficient. Even more freedom ex¬ 
ists for the free energy of interaction between the solute 
and the solvent, which needs to be modeled with some 
number of adjustable parameters. These parameters are 
to be tuned by matching the coarse-grained and micro¬ 
scopic descriptions. How to do this matching in practice 
remains an important open question. We proposed a spe¬ 
cific model with a single adjustable parameter here but 
it remains to be seen whether this model is appropriate 
on a case-by-case basis, and if not, to make adjustments 
to the equations by following the approach developed in 
this work. 

It is important to emphasize here that (11081) can be 
used to study colloidal suspensions of more than one 
colloidal particle, however, the description will only be 
accurate when the colloids are further than about one 
grid cell apart. This is because our modeling of the bare 
diffusion coefficient is based on equilibrium Green-Kubo 
expressions for a single particle. This will fail to give an 
accurate approximation when two particles come closer 
to each other than a grid spacing; at such short distances 
the hydrodynamic correlations among the diffusing parti¬ 
cles will not be captured accurately This is no different 
from prior work [H, [2^, [^, [37l - l4ll | where the hy¬ 


drodynamic interactions are only resolved up to at most 
the Stokeslet or Rotne-Prager level. Furthermore, when 
the nanocolloids come close to each other we expect that 
their direct interactions with the solvent molecules or 
with each other will be affected and other terms in (11081) 
will need to modified as well. Ultimately, as the density 
is increased there will be many nanoparticles per hydro- 
dynamic cell and in this case a coarse-grained theory of 
fluid mixtures should emerge. Such a theory could per¬ 
haps provide a bridge between macroscopic fluid mixture 
equations and dynamic density functional theories 
with hydrodynamic effects 

Finally, the present theory is isothermal as the energy 
density of the fluid is assumed to be a fast decaying vari¬ 
able as compared with mass and momentum variables. 
Of course, this precludes the study of thermal processes 
that arise in nanocolloidal suspensions in the presence 
of thermal gradients. The formulation of non-isothermal 
models is the subject of ongoing work. 
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Appendix A: Derivation of exact results 


In this appendix we obtain a number of exact results for the equilibrium probability distributions, the free energies, 
and the reversible drift term. 


1. Equilibrium distributions 


The equilibrium probability P{x) in Eq. (I4|) in the present level of description takes the form 

/o, g) = / dz-exp{-l3H{z)}S(R - qo) - p^(z))(5(g^ - gf,{z)) (Al) 

and the conditional expectation (• • • )^ in Eq. ([3|) takes the form 
Rpg ^ 


(. .. = 


p®q(R, 


/■ 1 . ^ 









25 


It is convenient to introduce the marginal equilibrium probability p), 

r \ ^ 

= J dz— exp{-PH{z)}S{R - qo)Y[S{pf, - p^{z)) 


(A3) 


and the corresponding conditional expectation (• ■ conditional on R, p, and not on g. Finally, we will also consider 
the equilibrium probability distribution of the density field, in the absence of nanoparticle, which is defined as 


N 


Psoiip) = I I[dq^dp,—^ 


exp 



l[Sip,-pf{z)) 


(A4) 


where is the normalization and the solvent mass density is introduced in Eq. (fT51) . The corresponding 

equilibrium expectation over solvent degrees of freedom conditional to the solvent density is denoted by (• • 


The three probabilities (lAll) . (IA3I) . (IA4I) are related to each other. By integrating the momentum variables in (lAll) . 
occurring in the kinetic energy of the Hamiltonian and in the Dirac delta functions, with (IE2I) in the Appendix lEl we 
have 


N 


M 


P®‘1(R,P,5) = — j y[dq*exp{-/3{7((7)}(5(R-qo)y[5(p^-p^(2;)) 


2=0 


exp|-§g^M^,)(z)g,.| 
(27r/,S)3W2 det M{z)^/‘^ 


(A5) 


where the mass matrix is defined in (j46l) . 


In a similar way, by integrating the atomic momenta in 
probability 


, gives the following form for the marginal equilibrium 


P'^'^(R,p) 



M 

{-/317(9)}<5(R-qo)n'^(PM 




(A6) 


where all momentum variables have been integrated out and Q is the normalization. We may write (IA5I) in the 
following form 


P‘>'^(R,P,5) = P«^(R,p) 


/ exp \ 

\ (27r//3)3W2 det M3/2 j 


(A7) 


that gives the relation between P®‘i(R, p, p) and P®'^ (R, p)- 


At the same time, the marginal P®‘i (R, p) in Eq. (IA3I) can be expressed in terms of the probability distribution of 
the solvent density Pso^(p) in the absence of nanoparticle. First, integrate momenta in (|A4p to get 


.N M 


(AS) 


where <5®°^ is the normalization. Then, Eq. ()A6p becomes 


N 


N 


M 


(R, p) = / n Q 1 (^'°‘(9)+E 1 f n - p 7 \^)) 


Z =1 


2=1 


= '^Psoi (p - moS{IL)) (|exp |-/3 J dr$“*(R - r) 


' sol 


p—mo(5(R.) 


sol 


exp{-/3$®^*(R)} 


(A9) 


where in the first equality we have integrated the Dirac delta function ^(R — qo) and in the second equality we have 
used the definition (|A8p . 
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2. Free energies 


The relationships (IA7I) and (IA9I) between the probabilities reflects into an exact expression for the free energy of 
the system. Let us introduce the free energy of the solvent the free energy of the fluid T(R.,p), and the CG 

Hamiltonian ^(R, p, g) (which is also a free energy) through the following expressions 

P;„'}(p)ocexp{-/3J-i(p)} 

(R, p) oc exp{-/3J^ (R, p) - 

P®'i(R, p, g) oc exp {-/3H(R, p, 5 )} (AlO) 


We now look at the relationships between these free energies (up to irrelevant constants). Because of (|A7| and 
(1X91) . the CG Hamiltonian has the form 


'H{R,p,g) 


—ksTln 


exp{-|g^M^^ig^} \ 
(27r//3)3“/2 detM3/2 j 


+ P(R, p) + $®’'*(R) 


(All) 


and the fluid free energy is 

P(R,p) =P"°‘ (p-mo(5(R)) +P“‘(R,p-mo5(R)) 
where the free energy of interaction between nanoparticle and solvent is 


J''"*(R,p) =-fcB'rin<J^exp<|-/3 / dr$'"*(R - r)h: 


:xsol 


P 

sol 


(A12) 


(A13) 


The exact result (IA12P that decomposes the free energy of the system into a solvent and an interaction part will be very 
useful for modelling. The fact that the free energy depends on p^ only through the combination p^°^ = p^i — mo<^At(R)) 
which is the mass density of the solvent, is a non-trivial result. 


3. The role of translation invariance on the free energy 


In this appendix we demonstrate the following exact identity involving derivatives of the free energy 

_ (R,,) . (A14) 

This identity is a direct consequence of translation invariance of the microscopic Hamiltonian. It is an important 
result because it gives an exact relationship between the derivatives of the free energy with respect to R and p^. This 
mathematical identity gives a strong condition on the modelling of the free energy. 

The proof is as follows. In the integrals over positions in (IA3I) perform the change of variables = q' — a which is 
a pure translation. The solvent potential is translation invariant and, therefore. 


N 


N 


(R, p) = y n ‘^'1* Q + E + a - q') + 

M / N \ 

X - ^m(5^(q' - a) 


(A15) 


2=1 


Note that the right hand side does not depend really on a. Take the derivative with respect to a and multiply by 
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ksT to obtain 


N 


N 


0 = / n 1 + E + a - + ‘i>""*(R) ) !> / rfrF“*(R + a - r)nr'(^) 


2=0 


2=1 


M 


AT 


(5 - mo(5^(R) - ^ m(5^(q' - a) 


N 


ksT J exp |-/3 ^ $“*(R + a - q') + $‘="‘(R) 


AT 


5 


M 


AT 


^E^y drV,5j.(r-a)p™'(z)]^,5 ( p^-mo(5^(R)-^m,5^(q'-a) 


(A16) 


Here, r'"*(R — r) is the force that a solvent particle located at r exerts on the nanoparticle located at R. Evaluate 
this expression at a = 0 and divide by P®'* (R, p) to obtain 

0 = J drr'"*(R-r)(hrE" + ^BT J dr^ Apeq (^17) 

Therefore, translation invariance implies the following exact result 

[ drF“‘(R - r) ^ ^soi^^^Rp ^ ^ 

J ^Pu ^Pv 

where we have used the notation (l44l) . This result is important because it relates the actual force on the nanoparticle 
due to the solvent with the derivatives of the free energy. 

Now, the gradient of the free energy F (R, p) introduced in Eq. (lAlOl) satisfies 


A^(R,rt + A4,»,R) = __^ 


N 


n exp |-/3 + E “ q*) “ 


N 


M 


N 


drF'^^(R-r)nr°H'^) Pfi -ruoS^iK) - 


N 


2 = 1 
N 


-ksT 


p-^(R,p)y 

r. M 


n dci,- exp -/3 u^°\q) + ^ $“‘(R - q.) - /3cI>®"*(R) 


2=1 

N 


(-moV(5,.(R)) - TOodp(R) - y^^mtS^{qi) 


(A19) 


i=l 


Hence we have the exact result 


4. The exact form of the reversible drift 


F(R) 


J drF'"*(R-r)(fi®°i)*^'’ 

dF dF 

-^(R,p)-moVd.(R)—(R,p) (A20) 


This result shows that the force on the nanoparticle due 
to the solvent is not simply the gradient of the free energy 
— ^ (R, p), but also depends on the chemical potential 
of the fluid near the particle. This is because the density 
includes the mass of the colloid in our formulation. The 
two exact result (IA18I) and (IA20I) combine to give the 
exact relation (|A14E 


We now consider the form of the drift given by (fTT|) 
The Poisson brackets entering the elements of the re- 
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versible matrix are computed as follows 

{r,r} = 0 
|r,p^| = 0 
{R,g^| = (5/,(ro) 

= 0 

r ^ ~ I — — 

Pti,Pu 2-j Qy.^ 


(A21) 


{PtJ.,S,^} = '^m,S„{ri)V6f,{r,) 

= J drpr{z)d„{r)V5^{r) 

= - '^miSf,{ri)VS^{r,) 
i 

= j drpr{z) 5 f,{T)V 5 y{Y) 

i 


= dr 


g“<5.(r)V^5^(r)-gf5^(r)V<5.(r) 

(A22) 


The conditional averages are 

({R,g.}) =5^(R) 

(A23) 


J 


where the double bracket notation is introduced in (l44ll . Therefore, the reversible part of the dynamics takes the 
form 


/ \ 


V (igp)"" J 

/ 

keT 


( 0 
0 


<5p(R) 


dn 


Jv V 



' dR ' 


OH 


dpu 


. 9 W . 


9gS 


9g? 


0 

.V“A 


y-V“5,(R)-^ 


^ ([g“<5._ |g/3j^v“,5,r^s 


(A24) 


where no approximations have been taken so far. 
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Appendix B: Approximate form for the reversible 
drift 


We will now use the approximations (l49)l ((52ll in order 
to compute all the different terms that appear in the 
exact equations (l43l) . Consider first the term 

= J dr (pr)"""® S^{r)VS,{r) (Bl) 

By using the linear for spiky approximation this be¬ 
comes 

/ dr{Mr)p.f^^S^{r)VS^{r) (B2) 

Note that the conditional expectation of the discrete den¬ 
sity field is just the conditioning value, this is = 

Per. Therefore, 


= J dr^,{v)p„5^V5^ = 




(B3) 


where we have used the definition of the interpolated den¬ 
sity field. By using the LFSA (1491) for the momentum 
field, the other required term becomes 




(B4) 


We see that, formally, the linear for spiky approximation 
approximate hatted functions with overlined functions, 
and allows to transform the double brackets |- • j^to 
simple space averages || • • • ||. 

Consider now the derivative of these terms that are 
required in (|A24I) 


^\\pS,V-6,\\ = 0 


(B5) 


This vanishes because the mass density and momentum 
density variables are independent. The next term is of 
the form 


d 

dpu 


d 




= = 0 (B6) 


where the term vanishes due to (lD7l) . Finally, we need 
to compute the following derivative 


d 




d 


= = 0 
where the terms vanish due to 


(B7) 


In summary, the form of the reversible dynamics under the linear for spiky approximation is 

/ \ / 0 0 S^iR) \ / m\ 


{Lp^i) 


R-pg 


V y 


0 


0 




dU 

dpu 


/ 


-ksT 


\-5,{R) -\p5,V<^5A \rdM5,\\-\g^5,V^5A ) \ J 


V -v“^^(R) y 


(B8) 


Let us now apply the linear for spiky approximation to 
the exact translation invariance identity (IA14I) . By mul¬ 
tiplying (IB3I) and ([_ 

p, by using (1291) . we obtain 


with the volume V)j, and sum over 


dpi 


IpV“<5.1^'>8 ||pV(5J| 


(B9) 


By using these approximations in the exact translation 
property (IA14|) we obtain the approximation 


dJ^ ^ 

as + “ “ 


(BIO) 


Appendix C: Approximate model for the CG 
Hamiltonian 


a. Modelling the kinetic part of the CG Hamiltonian 


The kinetic part of exact CG Hamiltonian in (IA11|) can 
be approximated under the LFSA (15^ in the form 


/ expj-fg^M ig,,| \ 

fesTln ( -^ ) 

\ (27r//3)3At/2 det M3/2 / 


Rp 


1 ttT-I 


ZkeT 


In det M 


(Cl) 


up to irrelevant constant terms. The order of magnitude 
of the term proportional to ksT can be estimate by as¬ 
suming a sufficiently smooth density field for which we 
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may approximate 

= (C2) 

y^i 

leading to a diagonal matrix. This approximation still 
satisfies the exact requirement m- The log det term of 
a diagonal matrix is simple 

j j. j ^ “iksT iksT 

—^— In det M = —^—tr In M. ~ ^^ In (C3) 

We observe that this term is not extensive, this is, does 
not scale as the number of particles per node. On the 
other hand, the kinetic energy 

2 

~ ^ ^ (C4) 

scales with the number of particles per node because, 
typically and p^ ~ giving 

(C5) 

which is an extensive quantity, proportional to the num¬ 
ber of particles per node. As we assume that the typical 
number of particles per node is large, we may neglect the 
term in det A4 in front of the kinetic energy term. 


From now on we will neglect this term and the CG 
Hamiltonian has the form 

H(R, P, g) = (R, p) + <i>-*(R) (C6) 

We will need the derivatives of the CG Hamiltonian that 
are given by 


dn _ dT 
M M dR 
dn 1 dM~Y' dT 

dp^ 2 ^ dpfj, dpf. 


dT. ^-1 


We now use the result 


dM^, 1 dM^,,,,,, 1 

- — 


dp^i 


dp^i 




(C7) 


(C8) 


become, finally 


m 

dR 

dn 

dPt, 

dn 

dgf. 




dR 


I — dT 

+ 71— 

dPi, 




MM 


(C9) 


b. Modelling the solvent part of the free energy 

The free energy of the solvent (p) is obtained from 
the first equation (lAlOl) from the probability (IA 8 I) . The 
explicit calculation of P^^^{p) is in general impossible due 
to the high dimensionality of the integrals in phase space. 
Therefore, we are forced to consider specific approximate 
models for this probability distribution. 

In accordance with the assumption that each hydrody¬ 
namic cell contains many fluid molecules, we will assume 
that the probability T’sq'|(p) is a Gaussian. The Gaussian 
probability has the form 

Psoiip) = exp (CIO) 

where N is the normalization, 6p^ = Pfj, — Peq are the 
fluctuations with respect to the homogeneous density peq, 
and the matrix of covariances is given by 

C^u = {hn^Pv)e^ = J dr J dr(5^(r)5i,(r') (5pr<5pr')eq 

(Cll) 

We estimate the form of this matrix as follows. We as¬ 
sume that the correlation of density fluctuations fluctua¬ 
tions decay in a length scale much smaller than the size 
of the cell and, therefore, the correlation can be approx¬ 
imated as proportional to the Dirac delta function, ac¬ 
cording to a standard result 

{SpJp^,)^^ = Y^S{r-r') (C12) 

where c is the isothermal speed of sound. The resulting 
free energy is quadratic in the density and will be termed 
Gaussian free energy. It has the explicit form 

= ^Sp^MPjp,. (C13) 

•^Peq 

This free energy function can be obtained from a lo¬ 
cal free energy functional of the form (square brackets 
denote a functional, while rounded parenthesis denote a 
function) 


where we have used the explicit form of the matrix in 
(155)) . With the discrete velocity (155)) the derivatives ()G7I) 


J^°'[p] = I drr\p{r)) 


(C14) 
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where /®°*(/o) is the thermodynamic free energy density 
of the solvent which, for the Gaussian model is 

fip) = :^ip- Peqf (C15) 

^Peq 

This functional is perhaps the simplest model famil¬ 
iar from Density Functional Theory. The model ne¬ 
glects molecular correlations, which is appropriate for 
the coarse description in which the hydrodynamic cells 
are much larger than molecular correlation lengths. The 
free energy (jC13ll is obtained from the functional (jC14ll - 
(IC15|) by using the interpolated field p(r) = in 

the functional, as advocated in Ref. |4^ . 

Once we have a free energy density, we may compute 
the pressure of the Gaussian model from the well-known 
thermodynamic relation 

jrsol 

p^°\p) =p^{p) - r\p) (C16) 

The pressure (IC16II that corresponds to (IG15I) is given by 
the quadratic equation of state (EOS) 

= (C17) 

^Peq 

Observe that for small deviations from equilibrium we 
obtain the expected linear EOS {p — Peq)- 


c. Modelling the interaction part of the free energy 

The interaction part of the free energy has the exact 
microscopic expression given in (I48II . We may obtain a 
simple model for this function if we consider the linear 
for spiky approximation (l49l) . Note that the microscopic 
potential energy of interaction between the nanoparticle 
and the solvent molecules can be expressed in terms of 
the microscopic solvent mass density as 

N 

y $“*(R - q,) = — / dr$“‘(R - r)y (z) (C18) 

^ m J 

1—1 

Within the linear for spiky approximation, we will ap¬ 
proximate the spiky field /5®°*(z) with a linear interpola¬ 
tion 

(C19) 

by using the approximation (IG19I) into (IG18I) we obtain 

^ 1 

y<l>“*(R-q,):^-ci>f(R)y(z) 

m ^ ^ 


where the nodal potential is defined according to 

$“(R) =y dr$“*(R-r)V’;x(r) (C21) 


We consider situations in which the nanoparticle is 
much smaller than the hydrodynamic cells and the range 
of the interaction potential $'"*(R — r) is also much 
smaller than the support of ^^(r). Therefore, we may 
approximate '0^(r) ~ i’fj.i'R) in Eq. (IG21I) . leading to 

4>”*(R) aV^^(R) (C22) 

where the constant a is the volume integral of the inter¬ 
action potential 

a= f dr$“*(r) = (C23) 

J Peq 

and we have introduced the “particle speed of sound” cq 
whose physical interpretation is that it gives the strength 
of the interaction of the nanoparticle with the solvent 
particles. Under these approximations, the microscopic 
potential of interaction between the nanoparticle and the 
solvent molecules is approximated by 

N 2 

y 4>”*(R - q,) (C24) 


Note that this approximation breaks translation invari¬ 
ance, because while the left hand side of (IG24I) is invariant 
under a translation of all the particles, the right hand side 
is not. In the approximation (IC24I) . the potential energy 
of the nanoparticle depends on the microscopic configu¬ 
ration z of the solvent particles only through the discrete 
solvent number densit y ff°^ {z). Therefore, by substitut¬ 
ing the approximation (IG24I) for the potential energy into 
the interaction part of the free energy (H51) we obtain the 
explicit model 

y) c. (C25) 

Peq 

By collecting (|C13|) and l|C25p the free energy (IA12I) be¬ 
comes 

p{R,p) = 

^Peq P^q 

(C26) 

In terms of the total mass density, the free energy is 

= -T^SPn^l^Jpu + 

^Peq Peq 

+ e(R) 


(C20) 


(C27) 
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where the last term is a density-independent term 


e(R) ^ ^ 


TT - Cn 




Peq _ 2 

The derivatives of the model (IC27I) are 


(C28) 


A 

9 




Peq 


T{B.,p) = —M^Jp,+ 

OPf, Peq 


mojcl - c^) 

Peq 


%{K) (C29) 


From Eq. (|A20I) the force on the nanoparticle due to the 
surrounding solvent is given by 


F(R) = (R, p) - moV5,(R)|^ (R, p) 

Peq 


(C30) 


This form of the force is consistent with the approxima¬ 
tion (|C24I) . 


We now consider the translation invariance (I66|) prop¬ 
erty of the free energy. This property is a strong guiding 
principle for the modelling of the free energy. In order to 
see if this important property is satisfied, we compute 


Wo (eg - c^) 


p^ (VV^^(R) + ||^^V(5.||V'.(R)) 


+ ^||pVp||+^(R) 

Peq V*—^ (yJrV 

=0 


mo(cg - c' 
Peq 


de 


(Vp(R) - IIAVpl)-f —(R) 


(C31) 


where we used that in a periodic domain ||pVp|| = 0 (as 
can be seen from integration by parts). Also, ||AVp|| is 
a compact notation for 

||AVp||=y'drA(R,r)Vp(r) (C32) 

and it depends on the position R of the nanoparticle. We 
observe that, in general, (|C31ll does not vanish. However, 
note that for sufficiently smooth density fields Eq (IS71) 
applies, and the first term is small. In particular, in 
the incompressible limit in which the density is constant, 
it vanishes identically. This strongly suggests that for 
modelling purposes, it is convenient to set e(R) = 0 and 
correct the free energy model developed so far in order 
to better respect translational invariance. 


solvent interacting with a single nanoparticle 

Peq 


T(R,p) = ^Sp^M^Jp. + 

■^Peq Peq 

(C33) 


Because this free energy gives the probability P(R, p), 
and we expect that for the case that the nanoparticle is 
identical to a tagged solvent particle this probability is 
Gaussian, we conclude that the limit of the nanoparticle 
becoming just another solvent particle is realized for co = 
c. 

The derivatives of the free energy model (IC33I) are 


A 

dR 


Peq 


^P(R, p) = —M^Jp, + ^ (C34) 

OPfj. Peq ^ Peq 

With these derivatives, we now compute the term (1551) 
entering the momentum equation 

JYi ^ ((A _ (A- 

=- \\S,75Vp\\ - ^ (^^(R)Vp(R) -(- ||p<5^VA|| 

Peq Peq 

= ~ (5^(R)Vp(R) - ||<5^AVp||) 


Peq 


(C35) 


where we have introduced the following total pressure 
equation of state 

F(r) = ^ {p(jf -Plq) + —^A(R,r)p(r) 

^Pec\ Peq 

(C36) 

Note that in the limit when the nanoparticle is just a 
tagged solvent particle we have cq = c and the last con¬ 
tribution to the pressure vanishes, giving simply the pres¬ 
sure of the Gaussian model. The last term in (jG35p is 
arguably small and will be neglected. Indeed, for smooth 
density fields 

||5^AVp|| ||5^A||Vp(R) = S^(R)Vp(R) (C37) 

and we have, finally 

= -||5;.VP|| (C38) 


1. Translation invariance and the barometric law 


In conclusion, in the present work we will use the fol¬ 
lowing model for the free energy of a fluid made of a 


In this appendix we examine the marginal probability 
P®'^(R) of finding the nanoparticle at position R. This 
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probability is, by definition, 

= J dpP^^{R,p) 

= y’dpexp{-/3(j-(R,p) + $""*(R))} (C39) 

Take its gradient 
B f 8 T 

—P'^'i(R) = -(3j dp—exp{-/3 (^(R,p) + <1>-‘(R))} 

n(F,ext 

- (C40) 

and use the approximate translation invariance of the free 
energy (pi 

8 a(f,ext 

^P-(R) + /3^P-(R) 

= pj dp||pVdJ|^exp{-^(j-(R,p) + <i>-‘(R))} 

= J dpexp{-^(^(R,p) + $-t(R))}A||pVd^|| 

= J dpexp{-pp}\\ip^V6^\\ = 0 (C41) 

Therefore, the (approximate) translation invariance 

property rigorously implies the well-known baromet¬ 

ric law 

P'=q(R) = 1 exp {-^$®’^*(R)} (C42) 

(qj 

where Q is the normalization factor. In the absence of 
an external field the probability to find the particle at a 
particular point R should be constant. 


Should the free energy model respect exactly the trans¬ 
lation invariance property (1661) . then the marginal distri¬ 
bution function would be rigorously given by the baro¬ 
metric law (IC42|) . However, the Gaussian model (IC33I) 
for the free energy satisfies (1661) only approximately, up to 
second order terms. As a consequence, the marginal dis¬ 
tribution P®'^(R) corresponding to the model (IC33I) does 
not give exactly the barometric law (IC42I) but rather 

P«q(R)= J dpexp{-/3P'(R,p)-/3$®’^*(R)} 

(X exp |-/3<i>-‘(R) + 

(C43) 

as can be seen by explicitly performing the Gaussian in¬ 
tegral. When = 0, the nanoparticle is not homoge¬ 
neously distributed in space but, rather, “sees” the un¬ 
derlying grid, unless it is a tagged fluid particle in which 
case Co = c. 


Appendix D: Derivatives of the basis functions 


In this work we assume periodic boundary conditions. 
Therefore any integration by parts give no surface terms. 
For example 

IIAVPil = -IIPVAII (Dl) 

for arbitrary functions A(r), B{r). 


We consider some identities that involve gradients of 
basis functions. For example, note that 

= 11(5^5. VV-. II = -\\^.V6^S4 

=- W-i^^Si^VS^W (D2) 

where in the second identity we have performed an inte¬ 
gration by parts. Therefore, 

||^/>^(5j.V5^|| =-2||(5^V^i,V(5^|| (D3) 

By multiplying both sides of this equation with and 

summing over /i, we obtain 

\\S.V^A=0 (D4) 

Another identity is obtained by introducing the following 
vector defined at each node p, 

(D5) 

Because of (lD4l) . this vector satisfies 

= 0 (D6) 

If the mesh of nodes is regular in such a way that for 
all nodes p we have a^ = ag, then the above equation 
implies ag = 0. Therefore, in regular grids we have the 
identities 

= 0 

||<5mV'.V,5,|| =0 (D7) 

It is expected that in non-regular meshes these quantities 
are also zero or very small. 
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Appendix E: Momentum integrals where the configuration dependent mass matrix is defined 

as 


In this appendix we quote the results for the following 
momentum integrals 


Iois)= I 


i=0 

M 


i=0 


N 






n^^o(27rw»fcBr)3/2 
(27rfcBr)3“/2 det M3/2 


exp<{ 


(El) 


N 

Mf,^{z) = ^m,(5^(qj)(5^(q*) (E3) 

i=0 

The above integrals are relatively easy to compute by us¬ 
ing the Fourier representation of the Dirac delta function. 


N f ^ 2 ) 

2=0 I 2=0 ^ J 

M / N \ 

X n ^ (X! - sa* j p* 

M \i=0 / 

= Io{g)m^S^{q,)M~^g^ (E2) 
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